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Synopsis : 


Concents : 


APPENDIX  I 


MILHY3  :  PROGRAMMING  DETAILS 


Laura  Baird 
University  of  Bristol 


May  I9S9 


Details  of  MILKY  computer  code.  The  subroutines 
are  documented  in  the  same  order  in  which  they 
appear  in  the  programme.  Functions  are 
documented  after  the  subroutine. 


Programme  changes  since  MILHY2 
MAIN 

Subroutines : 

HONDO 

STHYD 

RECHD 

CMPHYD 

SOILM 

HYDCON 

TWO 

GRAD 

SMCURV 

PRTHYD 

PUHYD 

HPLOT 

ADHYD 

SRC 

CMPRC 

STT 

CMPTT 

ROUTE 

RESVO 

ERROR 

SEDT 


Functions : 
GIT 
RMAX 
RMIN. 


BLOCK  DATA 


197 


PROGRAMME  CHANGES  SINCE  MILHY2 

I  have  made  major  changes  In  MILHY2  (coded  by  Dr  S  Howes,  University  of 
Bristol)  to  improve  the  predictive  capability  of  the  downstream  conveyance 
estimation,  specifically  under  out-of-bank  conditions.  ,  Therefore, 
turbulent  exchanges  of  momentum  between  segments  of  flow  in  a  cross-section 
and  multiple-routing  reaches  have  been  incorporated. 

However,  the  overall  structure  of  the  MILHY3  code  remains  unchanged  since 
MILHY2.  All  new  computations  occur  within  existing  subroutines,  some  of 
which  have  been  extensively  rewritten.  The  subroutines  which  have 
significantly  changed  are:- 

CMPRC 

ROUTE 

PRTHYD 

Minor  changes  have  occurred  in  subroutines :- 

CMpHYD 

ADHYD 

STHYD 

STT 

BLKDTA 

To  allow  for  the  nnew  computations,  BLKDTA  has  been  amended  to  allow  more 
data  to  be  read  from  'datal',  and  a  new  common  block,  '  COMMON’/ BL0CK2 '  has 
been  added  to  allow  this  new  data  to  be  transferred  between  subroutines. 


The  definition  of  some  arrays  in  ' COMMON/ BLOCK  1 '  have  been  changed.  These 
are 

A(  20 ,  70) 

Q( 20 ,  70) 

DEEP  (20,  70) 

C(  20 ,  6) 

■  .  ,  ,  ,  •  -  •  •  .  l-f 


end  area 
discharge 

elevation  of  water  surface 
absolute  elevation 
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The  array  size  of  A,  Q  and  DEEP  has  been  increased  to  allow  for  the 
generation  of  rating  curves  for  each  segment  of  flow  across  a  section, 
necessary  if  multiple  routing  reaches  are  invoked.  The  identification 
number  of  a  segment  is  firstly  the  cross-section  identification  number, 
followed  by  the  segment  number  (numbered  from  left  to  right  across  the 
section  looking  downstream). 

The  size  and  definition  of  the  array  C  (20,  6)  has  been  changed. 

It  was  previously  carried  in  'COMMON/BLOCK  1'  as  the  travel  time 
coefficient,  but  is  now  the  water  surface  elevation  plus  minimum 
cross-sectional  elevation.  The  travel  time  coefficient  is  now  carried  in 
COMMON/ BL0CK2  as  CC(20). 

ROIN  has  to  be  redefined  to  ensure  that  the  runoff  volume  is  correctly 
interpreted.  In  MILHY2  the  array  ROIN  was  the  hydrograph  volume  divided  by 
the  subcatchment  area,  hence  giving  the  depth  of  inundation  over  the 
catchment.  With  the  introduction  of  multiple  routing  reaches,  however,  it 
was  found  that  there  was  some  difficulty  in  carrying  the  correct  catchment 
area  due  to  the  repeated  ROUTE  command.  For  simplicity,  therefore, 
catchment  area  was  not  included  in  the  computation  of  ROIN  and  in  MILHY3 
ROIN  is  the  total  volume  of  discharge  (i.e.  the  area  under  the  hydrograph). 

Before  application  of  MILHY3  the  user  is  recommended  to  study  the  dataset 
details  for  data  1,  as  there  are  significant  changes  In  the  manner  in  which 
the  commands  need  to  be  entered  and  changes  in  the  order  of  the  data 
entries . 

The  programme  code  and  dataset  have  been  written  to  enable  the  user  to 
choose  which  of  the  conveyance  methods  are  most  suitable  for  the 
application.  O.uidelines  for  the  selection  of  these  methods  are  given  in 
the  final  report,  DAJA  . 
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SUBROUTINE  NAME: 
SYNOPSIS: 


COMMAND : 
INPUT: 


OUTPUT: 


VARIABLES  USED: 

Variables  held  in 

OCFS( 300 ,6) 

DATA( 311) 

CFS( 300) 

CTBLE( 50,11) 

RAIN( 300) 

ROIN(6) 

A(  20,70) 

Q( 20,70) 

DEEP(20, 70) 

ITBLE( 50,2) 

DP( 20) 

SCFS( 20) 

C(  20 , 6) 

ZALPHA( 20) 

IEND( 6) 

DA(  6 ) 

DISTC6) 

SEGN(6) 

DT(6) 

PEAK(6) 

ISG(6) 

NPU 

NHD 

NER 


MAIN 

Opens  files  for  input  and  output.  Initialises  certain 
variables.  Calls  appropriate  subroutine  according  to 
command . 


Two  data  files  called  'datal'  and  'data2'  must  exist. 

'datal'  -  contains  programme  controls  and  data 
'data2'  -  additional  information  for  the  infiltration 
algorithm 

See  appendix  X  and  II  for  details  of  these  two  data 
files 

Opens  the  output  file  'results'  to  which  the  details  of 
tie  simulation  are  to  be  sent. 


common  block  'BL0CK1': 
discharge 

data  input  for  each  command 
unit  hydrograph  discharge 
command  table 

cumulative  precipitation  at  equal  time  intervals 
volume  of  discharge  hydrograph 
end  area  (for  rating  curve) 
discharge  (for  rating  curve) 

elevation  of  water  surface  (for  rating  curve) 
integer  table 

flow  depth  for  previously  computed  travel  time 
flow  relationship 

discharge  for  previously  computed  travel  time 
flow  relationship 

absolute  stage  elevations  (for  rating  curve) 
alphanumeric  code  table 
number  of  points  in  the  hydrograph 
drainage  area 

segment  boundary  point  for  each  segment  of  a 
floodplain  and  chennel  cross-section 
Mannings  'n'  for  segment  of  a  floodplain  and 
channel  cross-section 

time  increment  for  precipitation  or  discharge 
data 

peak  discharge  of  hydrograph 

last  elevation  input  for  each  segment 

punch  code 

identification  number  of  hydrograph 
error 
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MAXNO 

NCOMM 

ICC 

NCODE 

TIME 

KCODE 

ICODE 


Variables  held 

PERQ(20, 70) 

TQ(20,6) 

CC( 20) 

LL(6) 

INRC 

LRC 


CONSTRAINTS: 


CALLED  BY: 

SUBROUTINES 

CALLED: 


maximum  number  of  data  entries  for  any  one  command 

number  of  command 

continuation  card 

command  number 

simulation  start  time 

measurement  unit  of  input  data 

if  (KC0DE=0)  then  imperial  units  else  metric 

required  measurement  unit  of  output 

if  (ICODE=0)  then  imperial  units  else  metric 

in  common  block.  '  BL0CK2  '  : 

percentage  discharge  for  segment  of  a  floodplain  and 
channel  cross-section  (for  rating  curve) 
total  discharge  for  cross-section  (for  rating  curve) 
travel  time  coefficient  for  previously  computed 
travel  time  relationship 

number  of  zero  discharge  values  for  segment  of  a 
floodplain  and  channel  cross-section 
identification  number  for  upstream  segment  rating 
curve 

identification  number  for  downstream  segment 
rating  curve 


The  most  important  constraints  in  this  program  involve 
the  limits  to  array  size  which  are  dimensioned  in 
COMMON.  For  example,  program  can  only  hold  6 
hydrographs  or  6  rating  curves  at  any  one  time.  17 
commands  which  are  defined  in  BLOCK  DATA  are  used  by 
MILHY3  (as  in  the  original  form  of  HYMO) .  The  legal 
MILHY  values  for  these  are  provided  in  HONDO  and 
appendix  I. 


HONDO 

STHYD 

RECHD 

CMPHYD 

PRTHYD 

PUHYD 

HPLOT 

ADHYD 

SRC 

CMPRC 

STT 

CMPTT 

ROUTE 

RESVO 

ERROR 

SEDT 


l _ £ 


FUNCTIONS 

NOTES 
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CALLED : 


There  are  changes  to  the  array  size  of  some  variables 
and  changes  to  the  definition  of  variables  in 
COMMON/ BL0CK1  from  MILHY2. 
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SUBROUTINE  NAME: 
SYNOPSIS: 

COMMAND: 

INPUT: 


OUTPUT: 


VARIABLES  USED: 


HONDO 

Reads  in  command  and  associated  data  from  file  'datal' 
and  by  comparison  to  the  legal  commands  contained  in 
CTBLE  (initialized  in  BLOCK  DATA),  it  determines  the 
command  number  (NCODE).  It  collects  up  the  variables 
from  the  variable  format  data  field. 


Data  is  read  in  from  file  'data!'. 

Command  must  be  located  in  the  first  20  columns  on  each 
line,  and  is  read  in  variable  ALPHA  (11)  (FORMAT 
2A1,9A2).  The  data  must  be  in  columns  21  to  80,  and  is 
read  into  variable  CHAR  (60)  (FORMAT  (60A1). 

Legal  commands  are: 

START 
STORE  HYD 
RECALL  HYD 
COMPUTE  HYD 
PRINT  HYD 
PUNCH  HYD 
PLOT  HYD 
ADD  HYD 

STORE  RATING  CURVE 
COMPUTE  RATING  CURVE 
STORE  TRAVEL  TIME 
COMPUTE  TRAVEL  TIME 
ROUTE 

ROUTE  RESERVOIR 
ERROR  ANALYSIS 
SEDIMENT  YIELD 
FINISH 

Additional  legal  entries  to  file  imlud": 

'*'  in  column  1  -  if  the  line  is  a  comment  line 

'*'  in  column  80  -  if  a  new  page  is  required  for  output 

Writes  out  the  command  and  associated  data  to  file 
'results',  and  returns  the  value  of  NCODE  to  MAIN  which 
is  then  used  to  select  the  next  subroutine  to  be 
called.  All  data  which  has  been  collected  for  the 
command  is  held  in  common,  in  the  array  DATA  (311). 


Variables  held  in  common  plus 


CHAR(60) 
ALPHA( 11) 
AUXA( 10) 
AUXB( 10) 


data  and  associated  text 
command 

array  used  to  collect  up 
array  used  to  collect  up 


data 
da  ta 
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CONSTRAINTS:  The  form  of  the  data  file  'datal'  must  be  strictly 

adhered  to.  HONDO  will  not  tolerate  spelling  mistakes. 

The  command  and  data  must  also  be  entered  into  the 
correct  columns.  The  data  must  be  in  the  order  which  is 
expected  by  HONDO  (these  details  are  provided  in 
appendix  I) . 

CALLED  BY:  MAIN 

SUBROUTINES 

CALLED: 

FUNCTIONS  GIT 

CALLED: 


NOTES 
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SUBROUTINE  NAME; 
SYNOPSIS: 

COMMAND : 

INPUT : 


OUTPUT: 


VARIABLES  USED: 


CONSTRAINTS: 

CALLED  BY: 

SUBROUTINES 

CALLED: 

FUNCTIONS 

CALLED: 


STHYD 

A  model  control  procedure.  Stores  the  coordinates  of  a 
measured  hydrograph  and  adds  a  constant  baseflow 
discharge  to  all  data  points. 


The  data  Input  for  this  command  has  been  read  into 
DATA(3ll)  by  HONDO  and  is  transferred  from  this  array 
into  the  following  variables  which  are  used  in  this 
subroutine : 


ID 

NHD 

DT(ID) 

DA(ID) 

BSF 

0CFS( 300 , ID) 

Stores  discharge  hydrograph  (time  and  associated 
discharge  values),  runoff  volume,  and  peak  discharge: 

OCFS(300,ID) 

ROIN(ID) 

PEAK(ID) 

Variables  held  in  common  plus 


ID 

NHD 

BSF 

DUMMY (300) 


storage  location  number 
hvdrograph  identification  number 
baseflow  discharge 
discharge  values  converted  to 
metric  units 


Only  6  hydrographs  at  any  one  time  can  be  stored  by  this 
program.  In  any  one  hydrograph,  a  maximum  of  300  points 
are  allowed. 


MAIN 


If  the  addition  of  baseflow  Is  not  required,  a  zero 
value  must  be  entered. 


NOTES: 
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SUBROUTINE  NAME: 
SYNOPSIS: 

COMMAND : 

INPUT: 


OUTPUT: 

VARIABLES  USED: 

CONSTRAINTS: 

CALLED  BY: 

SUBROUTINES 

CALLED: 

FUNCTIONS 

CALLED: 


RECHD 

Model  control  procedure. 

Recalls  previously  computed  and  punched  hydrographs. 
RECALL  HYD 

The  data  input  for  this  command  has  been  read  into 
DATA( 311)  by  HONDO  and  is  transferred  from  this  array 
into  the  following  variables  which  are  used  in  this 
subroutine : 

ID 

NHD 

DT(ID) 

DA(ID) 

PEAK(ID) 

ROIN(ID) 

IEND(ID) 

OCFS( 300 , ID) 

Stores  hydrograph  0CFS(300,ID) 

Variables  held  in  common  plus 

ID  storage  location  number 

NHD  Hydrograph  Identification  number 

Only  6  hydrographs  at  any  one  time  can  be  stored  by  this 
program.  In  any  one  hydrograph,  a  maximum  of  300  points 
are  allowed. 


MAIN 


NOTES: 
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SUBROUTINE  NAME :  CMPHYD 

SYNOPSIS:  Hydrological  procedure. 

Developes  unit  hydrograph  and  convolves  it  with 
incremental  runoff  to  produce  the  discharge  hydrograph. 
Runoff  is  derived  by  calling  subroutine  SOILM, 
utilizing  the  infiltration,  or  using  the  curve  number 
routine 

COMMAND :  COMPUTE  HYD 

INPUT:  Data  has  been  read  into  DATA(310)  by  HONDO  and  is 

transferred  from  this  array  into  the  following  variables 
which  are  used  in  this  subroutine: 

ID 

NHD 

DT(ID) 

DA(ID) 

CN 

HT 

XL 

RAIN( 300 ) 

OUTPUT:  Stores  the  calculated  discharge  hydrograph,  runoff 

volume,  and  peak  discharge 

0CFS( 300 , ID) 

ROI.V(ID) 

PEAK(ID) 

V.ARIABLES  USED:  Variables  held  in  common  plus 

ID  storage  location  number 

NHD  hydrograph  identification  number 

HT  difference  in  elevation 

XL  lengch  of  main  channel 

CONSTRAINTS:  A  maximum  of  6  hydrographs  can  be  stored. 

A  maximum  of  200  data  can  be  included  in  the 
precipitation  data. 

CALLED  BY:  MAIN 

SUBROUTINES  SOILM 

CALLED: 

FUNCTIONS 

CALLED: 

NOTES :  This  subroutine  has  been  updated  to  permit  the  user  to  select 

the  curve  number  routine  or  the  infiltration  algorithm. 
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SUBROUTINE  NAME: 
SYNOPSIS: 

COMMAND : 

INPUT: 


OUTPUT: 

VARIABLES  USED: 


SOILM 

Simulation  of  infiltration  and  hence  incremental  runoff 
associated  with  a  particular  storm  event,  and 
redistribution  of  soil  water  after  precipitation  ceases. 
Includes  a  stochastic  methodology  for  incorporating 
spatial  variability  of  soil  hydrological  properties. 


Certain  data  has  been  passed  from  CMPHYD  to  SOILM: 

DT(ID) 

IR 

CUMRAIN 


Remaining  variables  are  read  directly  into  SOILM  from 
data  file  'data2'.  Appendix  II  provides  the  details  of 
the  form  of  this  data  file  and  the  information  which  is 
required  by  SOILM. 


Provides  incremental  runoff  which  is  located  in 
DATA( 311)  and  which  is  passed  back  to  CMPHYD.  This 
runoff  is  at  the  same  time  interval  as  the  precipitation 
data  which  has  been  supplied  (DT(ID)). 


DT(ID) 

IR 

CUMRAIN( 251 ) 

TIME 

SIMDUR 

ALR 

AMR 

AF 

NLA 

NLB 


NL 

TC0M( 20) 
NSCOL 
IPCAREA 
SRI 

SR2 

SR  3 


Time  increment  for  precipitation  and 

hence  runoff  data 

number  of  rainfall  observations 

cumulative  rainfall  totals 

time  when  simulation  begins 

simulation  duration  (hours) 

rain  start  time  (hours) 

rain  stop  time  (hours) 

simulation  iteration  period  (seconds) 

number  of  cells  in  layer  1  in  soil 

column 

number  of  cells  in  layer  2  in  soil 


column 

number  of  cells  in  soil  column 

thickness  of  each  cell  (metres) 

number  of  soil  columns 

percent  are  occupied  by  soil  column 

soil  water  content  at  saturation, 

layer  1  in  soil  column 

soil  water  content  at  saturation, 

layer  2  in  soil  column 

soil  water  content  at  saturation, 

layer  3  in  soil  column 
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ASR1 

ASR2 

ASR3 

SSR1 

SSR2 

SSR3 

SATCON 

SATC0N2 

SATC0N3 

ASATCON 

ASATC0N2 

ASATC0N3 

SSATCONl 

SSATC0N2 

SSATC0N3 

DETCAP 

ADETCAP 

SDETCAP 

THETA( 20) 

ATHETA( 20) 

STHETA 

N’Q 

X(  20) 

Y(  20) 

X2(20) 

Y2(  20) 


Same  variable  definitions  as  the 

three  above,  but  variable  types 

are  DOUBLE  PRECISION  rather  than  REAL 

Standard  deviation  of  SRI 
Standard  deviation  of  SR2 
Standard  deviation  of  SR3 

saturated  hydraulic  conductivity 
(net'es  per  second)  layer  1 
saturated  hydraulic  conductivity 
(metres  per  second)  layer  2 
saturated  hydraulic  conductivity 
(metres  per  second)  layer  3 

Same  variable  definitions  as  the 
three  above,  but  variable  types  are  mean 
values  in  DOUBLE  PRECISION  rather  than 
REAL 

Standard  deviation  of  SATC0N1 
Standard  deviation  of  SATC0N2 
Standard  deviation  of  SATC0N3 

surface  detention  capacity  (metres) 

DOUBLE  PRECISION  surface  detention 
capacity 

Standard  deviation  of  detention 
capacity 

initial  soil  water  content  for  each 
cell  (cubic  metres  per  cubic  metres) 

DOUBLE  PRECISION  initial  soil  water 
content  (cubic  metres  per  cubic  metres) 
mean  value 

Standard  deviation  of  THETA(20) 

number  of  observations  on  soil  moisture 
characteristics  curve 
moisture  values  on  soil  moisture 
characteristic  curve  for  layer  1 
(cubic  metres  per  cubic  metres) 
suction  values  on  soil  moisture 
characteristic  curve  for  layer  1 
(metres) 

moisture  values  on  soil  moisture 
characteristic  curve  for  layer  2 
(cubic  metres  per  cubic  metres) 
suction  values  on  soil  moisture 
characteristic  curve  for  layer  2 
(metres) 
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X3(20)  moisture  values  on  soil  moisture 

characteristic  curve  for  layer  3 
(cubic  metres  per  cubic  metres) 

Y3(20)  suction  values  on  soil  moisture 

characteristic  curve  for  layer  3 
(metres) 

AX(20)  Same  variable  definitions  as  the 

AX2 (20)  X( 20) ,  X2( 20) ,  and  X3(20)  above, 

but  variable  types  are  mean  values 
AX3(20)  in  DOUBLE  PRECISION  rather  than  REAL 

SCURV1  Standard  deviation  of  soil  moisture 

characteristic  curve  for  layer  1 
SCURV2  Standard  deviation  of  soil  moisture 

characteristic  curve  for  layer  2 
SCURV3  Standard  deviation  of  soil  moistire 

characteristic  curve  for  layer  3 

EMAX  maximum  evaporation  during  the  day 

(metres  per  second) 

WT  write-out  time  interval  (hours) 

I0UT  determines  amount  of  output 

if  (I0UT=1)  total  output 
if  (IOUT=0)  shorter  output 

CONSTRAINTS:  A  maximum  of  10  soil  columns  for  any  one  subcatchment 

area  is  permitted. 

The  soil  moisture  characteristic  curve  can  be  defined  by 
up  to  a  maximum  of  20  points. 

The  soil  column  can  have  a  maximum  of  20  cells. 

The  Initial  soil  moisture  contents,  defined  for  each 
cell  at  the  start  of  simulation,  must  lie  within  the 
range  of  the  soil  moisture  characteristic  curve. 

CALLED  BY:  CMPHYD 

SUBROUTINES  HYDCON 

CALLED:  TWO 

CRAD 
SMCURV 

G05DDF(NAG  subroutine) 

FUNCTIONS  RMAX 

CALLED:  RMIN 


NOTES: 


SUBROUTINE  NAME:  HYDCON 


SYNOPSIS:  Calculates  hydraulic  conductivity  for  a  particular  layer 

in  the  soil  column  from  the  soil  moisture  characteristic 
curve,  using  the  Millington  and  Quirk  method. 

COMMAND: 

INPUT:  Variables  passed  from  SOILM: 

X(20) 

Y(20) 

SATCON 

SR 


OUTPUT: 

Unsaturated  hydraulic  conductivity  values  are  passed 

back  to 

SOILM  in  Z(20) . 

VARIABLES  USED: 

X(20) 

moisture  values  on  soil  moisture 
characteristic  curve  for  the  particular 
layer  (cubic  metres  per  cubic  metres) 

Y(20) 

suction  values  on  soil  moisture 
characteristic  curve  for  the  particular 
layer  (metres) 

SATCON 

saturated  hydraulic  conductivity  for 
the  particular  layer 

SR 

saturated  soil  moisture  content  for  the 
particular  layer 

Z(  20) 

unsaturated  hydraulic  conductivity 
values  corresponding  to  X(20)  above. 

CONSTRAINTS:  Maximum  points  on  the  soil  moisture  characteristic 

curve,  and  hence  the  hydraulic  function  is  20. 


CALLED  BY:  SOILM 

SUBROUTINES 

CALLED: 

FUNCTIONS 

CALLED: 

NOTES: 
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SUBROUTINE  NAME: 
SYNOPSIS: 

COMMAND : 

INPUT: 


OUTPUT: 


VARIABLES  USED: 


CONSTRAINTS: 


TWO 

Calculates  the  soil  water  pressure,  hydraulic  potential, 
and  hydraulic  conductivity  for  each  cell  in  the  soil 
column,  associated  with  a  particular  soil  water  content. 


Variables  passed  from  SOILM: 


NA 

N8 

G(20) 

GZ(20) 

Z(  20) 
X(20) 
Y<20) 
DEPTHS  20) 


Soil  water  pressure,  hydraulic  potential,  and  hydraulic 
conductivity  are  passed  back  to  SOILM 


SWP( 20) 
HP0T(20) 
C0ND( 20) 

NA 

NB 

THETA( 20) 
G(20) 


GZ( 20) 
X(20) 


Y(  20) 


2(20) 

DEPTH( 20) 

SWP(20) 

HPOTC20) 

COND(20) 


number  of  cells  in  layer  1 
number  of  cells  in  layer  2 
initial  soil  moisture  content  of 
each  cell 

gradient  of  soil  moisture  characteristic 
curve,  ie  grad  between  each  pair  of 
points 

gradient  of  hydraulic  function,  ie 
grad  between  each  pair  of  points 
moisture  values  on  soli  moisture 
characteristic  curve  for  the  particular 
layer  (cubic  metres  per  cubic  metres) 
suction  values  on  soil  moisture 
characteristic  curve  for  the  particular 
layer  (metres)  values 
unsaturated  hydraulic  conductivity 
values  corresponding  to  X(20)  above 
distance  from  surface  to  the  midpoint 
of  each  cell 

soil  water  pressure  of  each  cell 
hydraulic  potential  of  each  cell 
conductivity  of  each  cell 


A  maximum  of  20  ceils  in  the  soil  column  Is  permitted 


CALLED  BY: 


SOILM 
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SUBROUTINES 

COMMAND: 

FUNCTIONS 
COMMAND : 


fe 


NOTES: 


I 


1 
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SUBROUTINE  NAME: 
SYNOPSIS: 

INPUT: 

OUTPUT: 

VARIABLES  USED: 


CONSTRAINTS: 

CALLED  BY: 

SUBROUTINES 

COMMAND: 

FUNCTIONS 

COMMAND: 


GRAD 

Calculates  the  gradient  of  the  soil  moisture 
characteristic  curve,  and  hydraulic  conductivity 
function. 

Variables  passed  from  SOILM: 

X(  20) 

Y(20) 

Z(  20) 

Variables  containing  gradients  passed  back  to  SOILM. 

G(20) 

GZ(20) 

X(20)  moisture  values  on  soil  moisture 

characteristic  curve  for  the  particular 
layer  (cubic  metres  per  cubic  metres) 
Y(20)  suction  values  on  soil  moisture 

characteristic  curve  for  the  particular 
layer  (metres)  values 

Z(20)  unsaturated  hydraulic  conductivity 

values  corresponding  to  X(20)  above 
G(20)  gradient  of  soil  moisture  characteristic 

curve,  i.e.  gradient  between  each  pair 

of  points 

GZ(20)  gradient  of  hydraulic  function,  i.e. 

gradient  between  each  pair  of  points 

A  maximum  of  20  cells  in  the  soil  column  is  permitted 

SOILM 


NOTES: 
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SUBROUTINE  NAME: 
SYNOPSIS: 

INPUT: 


OUTPUT: 

VARIABLES  USED: 


CONSTRAINTS: 

CALLED  BY: 

SUBROUTINES 

CALLED: 

FUNCTIONS 

CALLED: 


SMCURV 


Generates  new  soil  moisture  characteristic  curve  based 
on  the  randomly  generated  moisture  values. 

Variables  passed  from  SOILM: 


AX(20) 

Y(20) 

SCURV 

SR 

NQ 

Coordinates  of  new  soil  moisture  characteristic  curve 
passed  back  to  SOILM: 


XNEW( 20) 
YNEW(20) 

AX(20) 

Y(  20) 

SCURV 

SR 

NQ 

XNEW( 20) 
YNEW( 20) 


values  of  soil  moisture  on  input  soil 
moisture  characteristic  curve  DOUBLE 
PRECISION  variable  type 

values  of  suction  on  input  soil  moisture 
characteristic  curve 
standard  deviation  of  soil  moisture 
characteristic  curve  in  DOUBLE  PRECISION 
saturated  soil  moisture  content 
number  of  coordinates  defining  soil 
moisture  characteristic  curve 
generated  soil  moisture  content  on  new 
soil  moisture  characteristic  curve 
generated  suction  values  on  new  soil 
moisture  characteristic  curve 


A  maximum  of  20  points  to  define  the  soil  moisture 
characteristic  curve 


SOILM 

G05DDF  (NAC  subroutine) 


RMIN 

RMAX 


NOTES: 
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SUBROUTINE  NAME; 
SYNOPSIS: 

COMMAND: 

INPUT: 


OUTPUT: 


VARIABLES  USED: 


PRTHYD 

Model  control  procedure. 

Prints  out  the  coordinates  of  a  hydrograph  and/or 
the  peak  value  and  runoff  volume. 

Converts  0CFS(300,ID)  to  a  stage  array,  S(300,ID) 
using  a  recalled  rating  curve. 

PRINT  HYD. 

The  data  input  for  this  commandn  has  been  read  into 
OCFS(300,ID)  by  HONDO  and  is  transferred  from  this  array 
into  the  following  variables  which  are  used  in  this 
subroutine: 

ID 

NPK 

IDR 

IN 


Details  of  the  hydrograph  are  held  in  common  and  are 
referenced  by  ID. 

Discharge,  DUMMY(300)  or  stage,  S300,ID)  hydrograph 
are  written  to  output  file  'results'. 

DUMMY( 300) 

S( 300 , ID) 

ROINI 

PEAK1 

PEAK5 

Variables  in  common  plus 


ID 

NPK 


IDR 


IN 


DUMMY ( 300) 

S( 300 , ID) 

PEAK1 

ROINI 

PEAK5 


storage  location  number 
form  of  output  required 
0  peak  and  volume  only 

1  discharge  hydrograph 

2  stage  hydrograph 

identification  number  of  rating  curve 
or  segment  to  be  used  for  conversion  to 
a  stage  hydrograph 
format  of  output 
0  five  columns  across  page 
1  single  column 

discharge  array  (converted  to  metric 
units  If  required) 

stage  array  (converted  to  metric  units 

if  required) 

peak  discharge 

volume  of  hydrograph 

peak  stage 
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CONSTRAINTS :  Maximum  of  300  points  define  the  hydrograph.  For 

conversion  to  stage  hydrograph,  rating  curve  must 
have  been  computed.  A  stage  hydrograph  cannot  be 
computed  if  multiple  routing  is  invoked. 

CALLED  BY:  MAIN 

SUBROUTINES 

CALLED: 

FUNCTIONS 

CALLED: 


NOTES: 

Conversion  to  a  stage  hydrograph  uses  a  previously  computed  rating 


217 


SUBROUTINE  NAME:  PUHYD 

SYNOPSIS:  Model  control  procedure. 

Punches  hydrograph  coordinates  on  cards  so  that  they 

can  be  recalled  by  RECHD 

COMMAND :  RECALL  HYD 

INPUT:  The  data  input  for  this  command  has  been  read  into 

DATA(311)  by  HONDO  and  is  transferred  from  this  array 
into  the  following  variables  which  are  used  in  this 
subroutine : 

ID 

Details  of  the  hydrograph  are  held  in  common  variables 
and  are  referenced  by  ID. 

OUTPUT:  Punched  version  of  0CFS(300,ID) 

VARIABLES  USED:  Variables  held  in  common  plus 

ID  storage  location  number 

DUMMY(300)  discharge  values  converted  to 

metric  units 

CONSTRAINTS:  Maximum  of  300  points  define  the  hydrograph. 

Maximum  of  6  hydrographs  in  store  at  any  one  time. 

CALLED  BY:  MAIN 

SUBROUTINES 

CALLED: 

FUNCTIONS 

CALLED: 


NOTES : 
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SUBROUTINE  NAME: 
SYNOPSIS: 

COMMAND: 

INPUT: 


OUTPUT: 

VARIABLES  USED: 

CONSTRAINTS: 

CALLED  BY: 

SUBROUTINES 

CALLED: 

FUNCTIONS 

CALLED: 


HPLOT 

Model  control  procedure. 

Plots  either  1  or  2  hydrographs  on  a  set  of  axis. 

PLOT  HYD. 

The  data  input  for  this  command  has  been  read  into 
DATA(311)  by  HONDO  and  is  transferred  from  this  array 
into  the  following  variables  which  are  used  in  this 
subroutine:- 

ID1 

ID2 

Details  of  the  2  hydrographs  are  held  in  common 
variables  and  are  references  by  ID1  and  ID2 

Discharge  plots  and  axis  are  written  to  output  file 
' results’ . 

CFS( 300) 

Variables  in  common  plus 

ID1 

ID2 

If  the  time  interval  of  the  two  hydrographs  to  be 
plotted  is  not  equal,  the  larger  increment  is  selected 
and  the  other  hydrograph  points  are  interpolated  at  this 
increment . 


MAIN 


NOTES: 
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SUBROUTINE  NAME: 
SYNOPSIS: 

COMMAND: 

INPUT: 


OUTPUT: 


VARIABLES  USED: 


CONSTRAINTS: 

CALLED  BY: 

SUBROUTINES 

CALLED: 

FUNCTIONS 

CALLED: 

NOTES: 


ADHYD 

Model  control  procedure 

Adds  together  the  coordinates  of  two  hydrographs 
ADD  HYD 

The  data  input  for  this  command  has  been  read  into 
DATA(311)  by  HONDO  and  is  transferred  from  this  array 
into  the  following  variables  which  are  used  in  this 
subroutine : 

ID 

NHD 

ID1 

ID2 

Details  of  the  2  hydrographs  are  held  in  common 
variables  and  are  referenced  by  ID1  and  ID2 

The  discharge  coordinates,  peak  discharge,  and  runoff 
volume  of  the  resultant  hydrograph: 

OCFS(300, ID) 

PEAK(ID) 

ROIN(ID) 


Variables  in  common  plus 


ID 

storage  location  number  for 
resultant  hydrograph 

NHD 

hydrograph  identification  number 
of  resultant  hydrograph 

ID1 ,  ID2 

storage  location  numbers  of  the 
two  hydrographs  to  be  added 

If  the  time  Interval  of  the  two  hydrographs  to  be  added 
Is  not  equal,  then  the  smaller  increment  is  selected  and 
the  other  hydrograph  points  are  interpolated  at  this 
increment . 


MAIN 


L 


A. 


SUBROUTINE  NAME: 
SYNOPSIS: 

COMMAND : 

INPUT: 


OUTPUT: 

VARIABLES  USED: 

CONSTRAINTS: 

CALLED  BY: 

SUBROUTINES 

CALLED: 

FUNCTIONS 

CALLED: 


SRC 

A  model  control  procedure 

Stores  a  rating  curve  in  form  of  elevation,  end  area, 
discharge  table 

STORE  RATING  CURVE 

The  data  input  for  this  command  has  been  read  into 
DATA(311)  by  HONDO  and  is  transferred  from  this  array 
into  the  followjng  variables  which  are  used  in  this 
subroutine : 

ID 

VS 

DEEP( 20 , ID) 

A(20, ID) 

Q(20, ID) 

Stores  the  rating  curve  in  variables  held  in  common: 

DEEP( 20 , ID) 

A( 20, ID) 

Q( 20, ID) 

Variables  held  in  common  plus 

ID  storage  location  number  of  rating 

curve 

VS  valley  cross  section  number 

Only  6  rating  curves  can  be  held  within  the  program 
at  any  one  time. 

Maximum  number  of  points  defining  rating  curve  are  20, 
MAIN 


NOTES: 
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SUBROUTINE  NAME: 
SYNOPSIS: 


COMMAND: 

INPUT: 


OUTPUT: 


VARIABLES  USED: 


CMPRC 

A  hydrological  procedure. 

Computes  rating  curve  for  valley  cross  section  using 
Mannings  equation.  If  turbulent  exchange  routines  are 
invoked  calculates  rating  curve  incorporating  momentum 
transfer  between  channel  and  floodplain  flows  during 
out-of-bank  conditions.  If  multiple  routing  reaches  are 
invoked  calculates  separate  rating  curves  for  each 
segment  of  the  cross-section  and  computes  the  percentage 
of  total  flow  which  would  occur  in  each  segment  at  the 
twenty  stage  computation  points. 

COMPUTE  RATING  CURVE 

The  data  input  for  this  command  has  been  read  into 
DATA(311)  by  HONDO  and  is  transferred  from  this  array 
into  the  following  variables  which  are  used  in  this 
subroutine : 

ID 

IT 

MR 

VS 

NSEG 

ELO 

EMAX 

SLOPE  1 

SL0PE2 

SEGN(NSEG) 

DIST(NSEG) 

DATA( 10:310) 

Stores  the  rating  curve  and  percentage  flow  In  each 
segment  in  variables  held  in  common 

A( 20, ID) 

Q( 20, ID) 

C(20, ID) 

DEEP( 20 , ID) 

PERQ(20, ID) 

TQ( 20, ID) 

Variables  heiu  in  common  plus 

ID  storage  location  number  for  rating  curve 

IT  turbulent  exchange  between  main  channel 

and  floodplains  invoked 

MR  multiply  routing  reaches  invoked 

VS  valley  section  identification  number 

NSEG  number  of  segments  in  valley  section 

ELO  lowest  elevation 
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EMAX  maximum  elevation 

SLOPE  1  channel  slope 

SL0PE2  flood  plain  slope 

DATA  (10:310)  alternate  distances  and  elevations 
(defining  cross  section) 

CONSTRAINTS:  Maximum  number  of  segments  in  a  cross  section  is  6. 

Maximum  number  of  points  in  a  rating  curve  is  20. 
Turbulent  exchange  and  multiple  routing  reaches  may  be 
invoked  independently  but  for  either  to  operate  there 
must  be  a  floodplain  segment  on  either  side  of  channel 
segments.  For  accuracy  the  user  is  recommended  to  have 
cross-sectional  positional  data  at  both  extremities  of 
flow  segments. 

CALLED  BY:  MAIN 

SUBROUTINES 

CALLED : 

FUNCTIONS 

CALLED: 


NOTES :  The  user  is  recommended  to  consult  MILHY3 :  datal  on 

the  form  of  dataset  'datal'  must  take  for  application 
of  multiple  routing  reaches. 
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SUBROUTINE  NAME: 
SYNOPSIS: 

COMMAND: 

INPUT: 


OUTPUT: 

VARIABLES  USED: 


CONSTRAINTS: 

CALLED  BY: 

SUBROUTINES 

CALLED: 

FUNCTIONS 

CALLED: 


STT 

Model  control  procedure. 

Stores  a  depth,  flow,  travel  time  table  (used  in  flood 
routing) . 

STORE  TRAVEL  TIME 

The  data  input  for  this  command  has  been  read  into 
DATA(311)  by  HONDO  and  is  transferred  from  this  array 
into  the  following  variables  which  are  used  in  this 
subroutine : 

ID 

REACH 

XL 

SLOPE 
MET1 
DP( 20) 

SCFS(20) 

CC(  20) 

Stores  travel  time  table  in  following  common  variables: 

DP( 20) 

SCFS(20) 

CC(20) 

Variables  held  in  common  plus 

ID  storage  location  number 

REACH  reach  identification  number 

XL  length  of  reach 

SLOPE  slope  of  reach 

A  maximum  of  20  points  are  allowed  to  define  a  travel 
time  table. 

MAIN 


NOTES: 
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SUBROUTINE  NAME: 
SYNOPSIS: 

COMMAND: 

INPUT: 


OUTPUT : 


VARIABLES  USED: 


CONSTRAINTS: 

CALLED  BY: 


CMPTT 

Hydrological  procedure. 

Compute  travel  time  table. 

COMPUTE  TRAVEL  TIME 

Data  has  been  read  into  DATA(3110)  by  HONDO  and  is 
transferred  from  this  array  into  the  following  variables 
which  are  used  in  this  subroutine: 

ID 

REACH 

NOVS 

XL 

SLOPE 

MR 

INRC 

LRC 

Stores  the  travel  time  table  in  following  common 
variables: 

DP(20) 

SCFS( 20) 

CC( 20) 


Variables  held  in  common  plus 


ID 

REACH 

NOVS 

XL 

SLOPE 

MR 

INRC 

LRC 


storage  location  number 

reach  identification  number 

number  of  valley  sections  in  the  reach 

length  of  reach 

Slope  of  reach 

multiple  routing  invoked 

upstream  segment  rating  curve 

identification  number 

downstream  segment  rating  curve 

identification  number 


A  maximum  of  20  points  are  allowed  to  define  a  travel 
time  table. 

A  maximum  of  6  valley  sections  are  permitted  in  a  reach, 
except  where  multiple  routing  reaches  are  invoked  where 
two  segment  section  must  be  identified. 


MAIN 


SUBROUTINES 

CALLED: 
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FUNCTIONS 

CALLED: 

NOTES: 


If  multiple  routing  reaches  are  invoked,  a  compute 
travel  time  table  and  route  command  must  be  entered  for 
each  segment  routing  reach. 
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SUBROUTINE  NAME :  ROUTE 

SYNOPSIS:  A  hydrological  procedure. 

Routes  a  hydrograph  through  a  reach  using  the  variable 
storage  coefficient  method. 

If  multiple  routing  reaches  are  invoked,  routes  a 
hydrograph  through  a  segment  routing  reach. 

Also  calculates  inflow  hydrograph  for  segment  reach 
using  PERQ{20)  from  rating  curve. 

COMMAND:  ROUTE 

INPUT:  The  data  input  for  this  command  has  been  read  into 

DATA(311)  by  HONDO  and  is  transferred  from  this  array 
into  the  following  variables  which  are  used  in  this 
subroutine : 

ID 

NHD 

IDH 

DT(ID) 

MR 

Details  of  the  hydrograph  to  be  routed  are  held  in 
common  variables  and  are  references  by  IDH. 

Details  of  the  inflow  segments  rating  curve  are  held  in 
common  variables  and  are  referenced  by: 

PERQ( 20) 

TQ(20) 

C(20, INRC) 

INRC. 

OUTPUT:  Stores  the  calculated  outflow  hydrograph,  its  peak 

discharge,  and  runoff  volume  in  common  variables: 

OCFS(300, ID) 

PEAK(ID) 

ROIN(ID) 

Proportional  discharge  for  each  time  increment  value  (if 
multiple  routing  reaches  are  invoked)  written  to  output 
file  ' results'  . 

D0CF(300, ID) 

VARIABLES  USED:  Variables  held  in  common  plus 

ID  storage  location  number  of  calculated 

outflow  hydrograph 

hydrograph  identification  number  of 
outflow  hydrograph 


NHD 
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CONSTRAINTS: 


CALLED  BY: 

SUBROUTINES 

CALLED: 

FUNCTIONS 

CALLED: 


IDH 

hydrograph 

hydrograph 

MR 

DOCFS 

P 


storage  location  number  of  inflow 
DT(ID)  iteration  period  of  outflow 

multiple  routing  invoked 
dummy  discharge  area  to  prevent 
overwriting  of  inflow  array 
percentage  of  inflow  in  multiple  routing 
reach  segment 


Discharges  included  in  the  inflow  hydrograph  must  be 
within  the  limits  of  the  travel  time  table,  otherwise 
there  is  no  way  to  define  the  travel  time  coefficient. 
If  the  solution  to  the  routing  equations  fails  to 
converge  after  10  Iterations,  convergence  is  forced. 

If  multiple  routing  reaches  are  invoked  the  inflow 
hydrograph  must  not  exceed  the  rating  curve  used  to 
compute  proportional  inflow  in  segment. 


MAIN 


NOTES :  If  multiple  routing  reaches  are  invoked,-  a  compute^ 

-  travel  time  table  and  route  command  must  be  entered  for 

each  segment  routing  reach.  Also  the  identification 
number  of  the  inflow  and  outflow  hydrographs  must  not  be 
the  same. 
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■  1  SUBROUTINE  NAME: 

SYNOPSIS: 

COMMAND: 

INPUT: 

< 

t 

. 

I 

> 

[ 

,  OUTPUT: 

VARIABLES  USED: 

CONSTRAINTS: 

,  CALLED  BY: 

l 

* 


RESVO 


A  hydrological  procedure. 

Routes  hydrograph  through  a  reservoir. 

ROUTE  RESERVOIR 

The  data  input  for  this  command  has  been  read  into 
DATA(311)  by  HONDO  and  is  transferred  from  this  array 
into  the  following  variables  which  a’-e  used  in  this 
subroutine : 

ID 

NHD 

IDH 

SCFS(20) 

STORE 

Details  of  the  inflow  hydrograph  are  held  in  common 
variables  and  are  referenced  by  ID: 

DT(ID) 

DA(ID) 

The  calculated  outdlow  hydrograph,  peak  discharge,  and 
runoff  volume  is  stored  in  common  variables: 

0CFS(300,ID) 

PEAK(ID) 

ROIN(ID) 

Variables  held  in  common  plus 


ID 

NHD 

IDH 

SCFS(20) 

STORE 


storage  location  number  of  calculated 
outflow  hydrograph 

hydrograph  identification  number  of 
outflow  hydrograph 
storage  location  number  of  Inflow 
hydrograph 

discharge  values  of  the  storage  discharge 
relationship  defined  for  the  reservoir 
storage  values  of  the  storage  discharge 
relationship  defined  from  the  reservoir 


The  discharge  of  the  inflow  hydrograph  must  be  within 
the  storage  discharge  relationship  defined  from  the 
reservoir.  A  maximum  of  20  points  are  allowed  to  define 
this  relationhip. 


MAIN 


1 _ g. 
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SUBROUTINES 

CALLED: 

FUNCTIONS 
CALLED : 


NOTES: 


SUBROUTINE  NAME:  ERROR 


SYNOPSIS: 

COMMAND : 
INPUT: 


OUTPUT : 


VARIABLES 


Model  control  procedure. 

Calculates  a  number  of  Indices  or  objective  functions 
which  detail  the  degree  of  fit  between  two  hydrographs. 

ERROR 

The  data  input  for  this  command  has  been  read  into 
DATA(311)  by  HONDO  and  is  transferred  from  this  array 
into  the  following  variables  which  are  used  in  this 
subroutine : 

ID1 

ID2 

Details  of  the  2  hydrographs  are  held  in  common 
variables  and  are  referenced  by  ID1  and  ID2 

The  values  of  HYMO's  original  error  statistics  plus  an 
additional  13  objective  functions  are  written  to  output 
file  ' results'  . 

ESDEV 
PCTER 
0F1 
OF  2 
0F3 
0F4 
OF  5 
0F6 
0F7 
OF8 
OF  9 
OF  10 
OF  1 1 

USED:  Variables  in  common  plus 


ID1 

storage  location  number  of  first 
hydrograph  (usually  assumed  to  be 
measured) 

ID2 

storage  location  number  of  second 
hydrograph  (usually  assumed  to  be 
calculated) 

ERR 

error  (measured  -  calculated  discharge) 

ESUEV 

error  standard  deviation 

PCTER 

percentage  peak  discharge  error 

0F1 

absolute  S'iro  of  errors 

0F2 

ordinary  least  squares 

0F3 

log  of  ordinary  least  squares 

OFi 

relative  sum  of  errors 
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CONSTRAINTS: 

CALLED  BY: 

SUBROUTINES 

CALLED: 

FUNCTIONS  CALLED: 


OF  5 

absolute  error 

difference 

0F6 

relative  error 

difference 

OF  7 

absolute  error 

divided  by  variance 

OF  8 

relative  error 

divided  by  variance 

OF9 

absolute  error 
by  variance 

difference  divided 

OF  10 

relative  error 
by  variance 

difference  divided 

0F1 1 

Pearsons  correlation  coefficient 

The  first  hydrograph  (ID1)  is  taken  to  be  the  measured. 
All  indices  are  printed  out  in  file  'results'  in  metric 
units . 

MAIN 


NOTES : 


SUBROUTINE  NAME:  SEPT 


SYNOPSIS: 

COMMAND: 

INPUT: 


OUTPUT: 

VARIABLES  USED: 


CONSTRAINTS: 

CALLED  BY: 

SUBROUTINES 

CALLED: 

FUNCTIONS 

CALLED: 


A  hydrological  procedure. 

Computes  Sediment  yield  for  a  f'eld  using  the 
Universal  soil  loss  equation. 

SEDIMENT  YIELD 

The  data  input  for  this  command  has  been  read  into 
DATA(311)  by  RONDO  and  is  transferred  from  this  array 
into  the  following  variables  which  are  used  in  this 
subroutine : 

ID 

SOIL 

CROP 

CP 

SL 

Details  of  the  hydrograph  are  held  in  common  variables 
and  are  referenced  by  ID: 

ROIN(ID) 

DA(ID) 

PEAK(ID) 

Writes  out  the  sediment  yield  to  the  output  file 
'results' : 

SED 

Variables  held  in  common  plus 

ID  storage  location  number  of  hvdrographs 

SOIL  soil  erodibllity  factor 

CROP  the  cropping  management  factor 

CP  the  erosion  control  practice  factor 

SL  the  slope  length  and  gradient  factor 


MAIN 


NOTES: 


T 
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FUNCTION  NAME: 

SYNOPSIS: 

INPUT: 


OUTPUT: 


GIT  (TCARD,  J,  JLAST,  SHIFT). 

Converts  alphabetic  array  to  floating  point  number. 


TCARD ( 10) 
J 

JLAST 
SHIFT 
A(  10) 

GIT 


VARIABLES  USED:  TCARD 
J 

JLAST 
SHIFT 
A(  10) 


CONSTRAINTS: 

CALLED  BY:  HONDO 


NOTES: 


FUNCTION  NAME: 


RMAX(X.NQ) 


SYNOPSIS: 

Returns 

the  maximum  element  in  a 

REAL  array 

INPUT: 

X(NQ) 

X  is  a  REAL  array 

size  NQ 

OUTPUT: 

RMAX 

VARIABLES  USED: 

X(NQ) 

RMAX 

CONSTRAINTS: 

CALLED  BY:  SOILM 

NOTES: 


FUNCTION  NAME: 

RMIN(X.NQ) 

SYNOPSIS: 

Returns  the 

minimum 

element  in  a  REAL  array. 

INPUT: 

X(NQ) 

X  is 

a  REAL  array  of  size  NQ 

OUTPUT: 

RMAX 

VARIABLES  USED: 

X(NQ) 

RMAX 

CONSTRAINTS: 

CALLED  BY: 

SOILM 

NOTES: 
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SUBROUTINE  NAME:  BLOCK  DATA 

SYNOPSIS:  Initializes  certain  variables.  These  variables  are  used 

to  determine  the  number  of  commands,  the  command,  the 
command  number,  the  maximum  number  of  data  entries  which 
are  associated  with  the  command,  and  the  data  which  is 
entered  in  the  variable  format  data  entry  area. 

INPUT: 

OUTPUT:  Initialized  arrays: 

ZALPHA( 20) 

CTBLE( 50,11) 

ITBLE( 50,2) 

NCOMM 

VARIABLES  USED:  ZALPHA(20)  alphanumeric  code  table  containing: 

1234567890  *., -(filled  with  blanks) 

CTBLE(50,11)  command  table  containing: 

START  (filled  with  blanks) 

STORE  HYD 
RECALL  HYD 
COMPUTE  HYD 
PRINT  HYD 
PUNCH  HYD 
PLOT  HYD 
ADD  HYD 

STORE  RATING  CURVE 
COMPUTE  RATING  CURVE 
STORE  TRAVEL  TIME 
COMPUTE  TRAVEL  TIME 
ROUTE 

ROUTE  RESERVOIR 
ERROR  ANALYSIS 
SEDIMENT  YIELD 
FINISH 

(filled  with  blanks) 

ITBLE(50,2)  integer  table  containing: 

I  4 
2310 
3310 
4310 

5  4 

6  1 

7  2 

8  4 
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9100 
10311 
11100 
12  8 

13  7 

14  25 

15  2 

16  5 

17  0 

NCOMM  number  of  commands  contains: 

7 


CONSTRAINTS: 

CALLED  BY: 
SUBROUTINES  CALLED: 
FUNCTIONS  CALLED: 


NOTES: 
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APPENDIX  2  datal 


The  commands  permitted  by  MILHY3  are  the  same  as  those  for  the  original  MILHY 
They  perform  two  functions:  model  control  and  hydrological  procedures.  The 
legal  commands  are: 

Model  control  procedures: 

START 

STORE  HYD 

RECALL  HYD 

STORE  RATING  CURVE 

STORE  TRAVEL  TIME 

ADD  HYD 

PRINT  HYD 

PLOT  HYD 

PUNCH  HYD 

ERROR  ANALYSIS 

FINISH 


Hydrological  procedures: 

COMPUTE  HYD 
COMPUTE  RATING  CURVE 
COMPUTE  TRAVEL  TIME 
ROUTE 

ROUTE  RESERVOIR 
SEDIMENT  YIELD 


An  *  in  column  1  means  that  the  line  is  a  comment. 

An  *  in  column  80  means  skip  to  a  new  page  before  writing  to  file. 


The  command  must  be  in  columns  1-20,  and  the  data  in  21-80. 


The  variables  in  the  data  area  MUST  be  in  the  correct  order,  but  can  be 
surrounded  by  text  (i.e.  units  and  variable  labels). 


The  following  two  pages  list  the  variables  which  are  required  for  each 
command,  and  the  order  in  which  they  are  required. 


Data  requirements  of  HYM03 
(data  to  be  located  In  'data  1') 


HYDROLOGICAL  PROCEDURE 


Variable 


COMPUTE  HYD 


Storage  location  number  for  hydrograph 

Hydrograph  identification  number 

Time  increment  for  rainfall  data  (hours) 

Watershed  area  (sqmi/km  j 

Curve  number  (enter  zero  if  not  invoked) 

Watershed  height,  maximum  difference  (ft/m) 

Main  stream  length  (mi/km) 

Rainfall,  cumulative  totals  (inches/mm) 


COMPUTE  RATING  CURVE 


Storage  location  number  for  rating  curve 
Turbulent  exchange  of  momentum  between 
segments  (not  invoked  enter  zero) 
(invoked,  enter  1-4  depending  on  method) 
Multiple  routing  reaches 

(not  invoked  enter  zero) 

(invoked  enter  one) 

Valley  section  location  number 

Number  of  segments  in  channel  (max.  of  6) 

Minimum  elevation  (ft/m) 

Maximum  elevation  (ft/m) 

Channel  slope 

Floodplain  slope 

Manning  'n'  for  each  segment 

(negative  value  for  channel  segments) 
Segment  boundary  points  (horizontal  distance) 
(ft/m) 

Cross-section  co-ordinates  (distance  then 
elevation)  (ft/m) 


used  in  subroutine 


ID 

NHD 

DT(ID) 

DA(ID) 

CN 

HT 

XL 

RAIN(300) 


ID 

IT 


MR 


VS 

NSEG 

ELO 

EMAX 

SLOPE  1 
SL0PE2 
SEGN(NSEG) 

DIST(NSEG) 

DATA( 12,311) 


/Compute  Travel  Time 
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COMPUTE  TRAVEL  TIME 


Storage  location  of  travel  time  table  ID 

Reach  identification  number  REACH 

Number  of  valley  sections  in  reach  NOVS 

Reach  length  (ft/m)  XL 

Slope  (average  for  flow  segments)  SLOPE 

Multiple  routing  reaches  MR 

(not  invoked  do  not  enter) 

(invoked  enter  one) 

Inflow  rating  curve  identification  INRC 

Outflow  rating  curve  identification  LRC 


(first  digit  is  storage  location  of  the 
rating  curve,  the  second  digit  is  the 
segment  number) 

N.B.  If  multiple  routing  reaches  not  invoked 
do  not  enter  values  for  INRC  and  LRC 


ROUTE 


Storge  location  number  of  outflow  hydrograph  ID 

Hydrograph  identification  number  of  outflow  NHD 

hydrograph 

Storage  location  number  of  inflow  hydrograph  IDH 

Time  increment  (hrs)  DT(ID) 

Multiple  routing  reaches  MR 

(not  invoked,  do  not  enter) 

(invoked  enter  one) 


ROUTE  RESERVOIR 


Storage  location  number  of  outflow  hydrograph  ID 
Hydrograph  Identification  number  of  outflow  NHD 

hydrograph 

Storage  location  number  of  inflow  hydrograph  IDH 
Reservoir  outflow  storage  relation  (max  20  DT(ID) 

points) 


SEDIMENT  YIELD 

Storage  location  of  number  of  hydrograph  ID 

Soil,  crop,  conservation  and  gradient  factors  SOIL,  CROP,  CP,  SL 


/Model  Control  Procedures 
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MODEL  CONTROL  PROCEDURES 


Start 


Start  time  (hours) 

Punch  code 

(not  invoked  enter  zero) 
(invoked  enter  one) 

Data  input 

imperial  enter  zero 
metric  enter  one 
Data  output 

Imperial  enter  zero 
metric  enter  one 


STORE  HYD 


Storage  location  number  for  hydrograph 
Hydrograph  identification  number 
Time  increment  for  discharge  data  (hrs) 
Watershed  area  (sq.mi/kra  ) 

Baseflow  (added  to  discharge)  (cfs/m_s 
Discharge  (300  points  max.)  (cfs/nrs  i) 


) 


RECALL  HYD 


Storage  location  number  for  hydrograph 
Hydrograph  identification  number 
Time  increment  for  discharge  data  (hours) 
Watershed  area  ( sq .mi, /km4) 

Peak  discharge  (cfs/m  s  A) 

Runoff  volume  (cf/ni  ) 

Number  of  points  in  hydrograph 
Discharge  (cfs/mJs  A) 


STORAGE  RATING  CURVE 


Storage  location  number  for  rating  curve 
Valley  section  number 
Rating  curve  points 
elevation  (ft/ml 
end  area  (ftz/ni  1 
flow  rate  (cfs/nrs  l) 


A 


TIME 

NPU 


KCODE 


ICODE 


ID 

NHD 

DT(ID) 

DA(ID) 

BSF 

OCFS(300 , ID) 


ID 

NHD 

DT(ID) 

DA(ID) 
PEAK(ID) 
ROIN(ID) 
IEND(ID) 
0CFS( 300 , ID) 


ID 

VS 

DEEP( 20 , ID) 
A(20 , ID) 
Q(20 , ID) 


/Storage  Travel  Time 
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STORAGE  TRAVEL  TIME 


Storage  location  number  for  travel  time  table  ID 

Reach  identification  number  NHD 

Length  of  reach  (ft/m)  XL 

Slope  either  channel  or  flood  plain  or  SLOPE 

weighted  average  of  the  two 

Depth  (ft/m)  .  DP(ID) 

Discharge  (cfs/nTs  X)  SCFS(20) 

Storage  coefficient  C(20) 


ADD  HYP 

Storage  location  numfe’.  for  resultant 
hydrograph 

Hydrograph  identification  number  of 
resultant 

Storage  location  of  two  hydrographs 
to  be  added 


PRINT  HYP 

Storage  location  number  of  hydrographs  ID 

Specification  of  type  of  output  NPK 

3  —1 

0  peak  and  volume  only  (cfs/m.s  ) 

1  discharge  hydrograph  (cfs/nrs  j 

2  stage  hydrograph  (ft/m) 

Rating  curve  identification  for  conversion  of  1DR 
hydrograph 

PLOT  HYP 

Storage  location  number  of  the  1  or  2  ID1,  ID2 

hydrographs  to  be  plotted 

PUNCH  HYP 

Storage  location  number  of  hydrograph  ID 

ERROR  ANALYSIS 

Storage  location  numbers  of  2  ID1 ,  ID2 

hydrographs  to  be  compared 

FINISH 


ID 

NHD 

ID1,  ID2 


No  information  required 
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Setting  up  datal:  changes  from  MILHY2 

Listed  in  the  previous  few  pages  are  the  legal  commands  and  variables  for 
dataset  'datal'.  Whilst  the  commands  remain  unchanged  the  attention  of  the 
user  is  drawn  to  order  of  the  variables  and  the  introduction  of  new  variables. 
These  new  variables  are  all  model  control  variables  allowing  the  user  to 
select  the  method  of  incorporating  turbulent  exchange  between  cross-sectional 
flow  segments  and  to  introduce  multiple-routing  reaches. 

If  the  multiple  routing  reaches  or  the  turbulent  exchange  algorithms  are 
invoked  along  a  particular  reach  then  any  channel  cross-sectional  segment  in 
that  reach  must  have  floodplain  segments  on  either  side.  In  addition,  if 
multiple  routing  is  invoked  there  may  be  only  two  cross-sections  along  the 
reach.  After  the  rating  curves  have  been  computed  for  these  two  cross 
sections  each  segment  of  flow  must  be  separately  routed.  This  means  there 
must  be  the  same  number  of  segments  in  both  the  upstream  and  downstream 
cross-sections.  For  each  segment  a  separate  'TRAVEL  TIME'  and  'ROUTE1  command 
are  entered  and  the  outflow  hydrographs  from  segment  are  added  to  give  the 
total  discharge  across  the  cross  section.  Care  must  be  taken  to  ensure 
storage  location  number  (ID's)  of  inflow  hydrographs  are  not  reused  in  the 
application  of  multiple-routing  reaches.  The  variables  INRC  and  LRC  in  the 
'TRAVEL  TIME'  command  identify  the  inflow  and  outflow  segments  for  one  reach 
of  a  multiple-routing  reach.  These  identification  numbers  consist  of  the 
cross-section  ID  number  followed  by  the  segment  number,  numbered  left  to  right 
looking  downstream.  In  the  ROUTE  subroutine  the  inflow  hydrograph  is 
distributed  between  the  segments  of  the  upstream  cross-section  using  the 
rating  curve  identified  by  INRC.  It  is  important  therefore  that  the  upstream 
cross-section  is  positioned  at  or  near  the  upstream  extremity  of  the  reach. 

An  example  'datal'  dataset  is  given  for  a  single  reach  where  both  multiple 
routing  and  turbulent  exchange  algorithms  are  Invoked. 


MILHY3  -  a  mathematical  flood  forecasting  model  for 
ungauged  catchments 


(MILHY2)  including  two-stage  channel  modelling. 

With  improved  out-of-bank  flood  modelling  incorporating 
TURBULENT  EXCHANGE  between  in  and  out  of  bank  flows  and 
MULTIPLE  ROUTING  REACHES  -allowing  separate  pathways  for 
channel  and  floodplain  flows. 


C  Coded  by: 
C 


University  of  Bristol 


Altered  subrout -es 


The  structure  of  the  code  remains  unaltered. 

All  additional  computations  occur  within  existing 
subroutines . 

HOWEVER,  there  are  significant  changes  in  the  manner  in 
which  the  dataset  DATA1  must  be  set  out  to  facilitate 
utilisation  of  the  new  capabilities. 

All  computations  (except  in  the  infiltration  algorithm) 
are  carried  out  in  imperial  units,  irrespective  of  KCODE 
and  ICODE.  C0*i10N/ BLOCK  1  and  CC*MDN/BL0CK2  use  imperial. 


C0M10N / BL0CK1 /  OCESC 300 , 6) , DATA(311) ,CrS( 300 ) , CTBLE( 50,11), 
&RAIN(300 ) ,ROIN(6) , 

&A(20,70) ,Q(20,70) ,DEEP(20, 70), ITBLE( 50,2), DP( 20), SCFSC20) ,C(20, 6) , 
&ZALPHA( 20 ) , I END ( 6 ) ,DA(6) ,DIST(6> ,SEGN(6> ,DT(6> ,PEAK(6) , ISG(6) , 
&NPU ,  NHD  ,  NER ,  MAXNO ,  NCOW ,  ICC ,  NCODE ,  TIME ,  KCODE ,  ICODE 
COf'MDN / BL0CK2 /  PERQ(20 , 70 ) .  TQ(20 , 6)  ,CC(20 )  .  LL(6 ) ,  INRC  , LRC 


C  Definition  of  variables  in  common  1 


Hydrograph  discharge 
Data  associated  with  each  conmand 
Unit  hydrograph  discharge 
Command  table 

Cumulative  precipitation  values 
Volume  of  discharge  hydrograph 

note:  this  variable  is  no  longer  divided  by  area 
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p 

i 


I 


C  A  End  area 

C  Q  Flow  rate  for  rating  curve 

C  DEEP  Elevation  of  water  surface  (for  rating  curve) 

C  ITBLE  Integer  table 

C  DP  Flow  depth  for  previously  computed  travel  time  flo*  relationship 

C  SCFS  Discharge  for  previously  computed  travel  time  flow  relationship 
C  C  Absolute  stage  elevations  computed  in  rating  curve 

C  ZALFA  Alphnumeric  code  table 
C  IEND  Number  of  points  1.'.  the  hydrograph 
C  DA  Drainage  area 

C  DIST  Segment  boundary  point  for  each  segment  of  a  cross  section 
C  SEGN  Mannings  'n'  for  each  segment  of  a  cross  section 
C  DT  Time  increment  for  rainfall  or  discharge 

C  PEAX  Peak  discharge  for  hydrograph 

C  ISG  Last  elevation  input  in  each  segment  position 
C  NPU  Punch  cod# 

C  NHD  Hydrograph  identification  number 

C  NER  Error  number 

C  MAXNO  Maximum  number  of  data  entires  to  be  expected  for  any  consnand 
C  NCCm  Number  of  coomands 

C  ICC  Continuation  line 

C  NCODE  Number  of  corrmand 

C  TIME  Start  time  of  simulation 

C  KCODE  Measurement  unit  of  input 

C  0  -  imperial 

C  not  0  -  metric 

C  ICODE  Measurement  unit  of  output 
C  0  -  imperial 

C  not  0  -  metric 

C  Variables  cocimon  2 

C  PERQ  Percentage  discharge  (rating  curve  computation) 

C  TQ  Total  discharge  (rating  curve  computation) 

C  CC  Travel  time  coefficient 

C  LL  Number  of  zero  discharge  values  in  rating  curve  segment 

C  INRC  Inflow  rating  curve  (multiple  routing) 

C  LRC  Outflow  rating  curve  (multiple  routing) 


1 

2 

3 


OPEN  (1, STATUS-' old1 .FILE-’ datal ’ ) 

OPENU5,FILE-’data2'  .STATUS-’ old*  ) 

OPEN (6. STATUS- 'new1 . FILE-' results ’ > 

NCODE- 0 
NPU-0 
ICC-0 
NER-0 

CALL  HONDO 
IF  (NER)  2,2.19 

GO  TO  (3.4. 5,6,7. 8, 9, 10. 11.12,13.14.15.16,17,18.19).  NCODE 
TIME-DATA(l) 

NPU-DATA(2) 

KCODE-DATA ( 3 ) 

IC0DE«DATA(4) 

GO  TO  1 
CALL  STHYD 
GO  TO  1 
CALL  RECHD 
GO  TO  l 
CALL  CMPHYO 
GO  TO  1 


6 
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7 

CALL  PRTHYD 

GO  TO  1 

8 

CALL  PUHYD 

GO  TO  1 

9 

CALL  HPLOT 

GO  TO  1 

10 

CALL  ADHYD 

GO  TO  1 

11 

CALL  SRC 

GO  TO  1 

12 

CALL  CMFRC 

GO  TO  1 

13 

CALL  STT 

GO  TO  1 

14 

CALL  CMPTT 

GO  TO  1 

15 

CALL  ROUTE 

GO  TO  1 

16 

CALL  RESVO 

GO  TO  1 

17 

CALL  ERROR 

GO  TO  1 

10 

CALL  SEDT 

GO  TO  1 

19 

STOP 

END 

SUBROUTINE  HONDO 


C  This  subroutine  reads  in  the  data  from  'datal',  searches  an  alphanumeric 
C  code  table  to  determine  the  NCODE  of  the  required  operation,  and  collects 
C  variables  from  the  freefloating  data  field. 

C  The  coomand  table  (CTBLE),  integer  table  (ITBLE),  number  of  coomands 
C  (NCOffl)  and  alphanumeric  array  (ZALPHA)  are  initialized  in  BLOCK  DATA 
C  Located  at  the  end  of  this  listing. 


COWON / BLOCK  1  /  OCFS(  300 . 6 ) ,  DATA(  3 1 1 ) .  CFS  ( 300 ) ,  CTBLE  ('0.11). 
&RAIN(300) ,ROrN(6) , 

&A( 20 , 70) ,Q( 20 , 70) , DEEP <  20 , 70 ) , ITBLE ( 50 , 2 ) , DP(2Q ) , SCFS ( 20 ) , C ( 20 . 6 ) , 
&ZALPHA(20 ) , IEND( 6 ),DA(6),DIST(6), SEGN( 6 ) , DT (6) , PEAK (6),ISG(6), 

&NPU . NHD, NER, MAXNO . NCOW .  ICC . NCODE ,  UME . KCODE .  ICODE 
COWON / BLOCK2 /  PERQ<20 . 70) . TQ<20 , 6) ,CC(20 ) ,LL(6) , INRC.LRC 

DIMENSION  CHAR( 60 ) ,  ALPHA ( 11 ) , AUXA( 10 ) , AUXB( 10 ) 

IF  (ICC)  1,1,3 
C  READ  IN  DATA  CARD 

1  READ  (1,42)  ( ALPHA ( I ) . 1*1 . 11 ) , (CHAR( I ) , I“1 . 60 ) 

C  IF  FIRST  CHARACTER  IS  BLANK  THE  CARD  IS  A  CONTINUATION  OF 
C  PREVIOUS  CARD . 

IF  (ALPHA(l)-EALPHA(ll))  2,9,2 

2  IF  (ICC)  3.3,40 

C  ASTERISK  IN  COL.  80  MEANS  SKIP  TO  NEW  PAGE  BEFORE  PRINTING  CARD 

3  IF  ( CHAR  ( 60  )  -  2ALPHA  (ID)  4,5.4 


4  WRITE  ($,43) 

5  WRITE  (6,44)  (ALPHA(I).I-l, 11), (CHAR(I), 1-1,60) 

C  IF  FIRST  CHARACTER  IS  A  *  THE  PREVIOUS  CARD  WAS  A  CCWEHT  CARD 
IF  ( ALPHA ( 1)-ZALPHA( 12) )  10,6,10 
C  IF  PUNCH  CODE  POSITIVE,  COHENT  CARDS  ARE  PUNCHED. 

6  IF  (NPU)  8,8,7 

7  WRITE  (7,45)  (ALPHA(I) ,1-1,11), (CHAR(I) ,1-1,60) 

8  ICC-0 
GO  TO  1 

9  WRITE  (6,44)  (ALPHA(I),I-1, 11). (CHAR(I). 1-1,60) 

GO  TO  24 

C  SEARCH  FIRST  TVtt  ALPHAMERIC  CHARACTERS  TO  SEE  IF  THEY  ARE  NUMBERS 

10  ICC-1 

DO  12  1-1,10 

IF  (ALPHA(l)-ZALPHA(I))  11.15.11 

11  IF  (ALPHA(2)-ZALPHA(I))  12,15,12 

12  CONTINUE 

C  STATEMENT  NUMBER  7  IS  BRANCHED  TO  IF  NUMBERS  ARE  PRESENT 
C  IF  NOT  NUMBER  SEARCH  COMMAND  TABLE  FOR  MATCH 

C  CALL  FIRST  10  VALUES  FROM  PERMANENT  DATA  STORAGE 

DO  14  1-1 ,  NCOti 
DO  13  J-1,11 

IF  (CTBLE(I, J)-ALPHA(J))  14.13,14 
C  SN  10-PART  MATCH 

13  CONTINUE 

C  IF  THIS  LOOP  IS  COMPLETED  WE  HAVE  COMPLETE  MATCH-  CALL  NCODE 

c  and  max  number  and  exit  loop 

NCODE-ITBLE (1,1) 

MAXNOITBLE(I,2) 

GO  TO  21 

14  CONTINUE 

C  IF  MAJOR  LOOPS  FINISHED  WITHOUT  A  MATCH  WRITE  ERROR  MESSAGE 
C  AND  SET  NER  -  1 

NER-1 

WRITE  (6,46) 

RETURN 

C  CONVERT  DIGIT  INPUT  CODE  FROM  ALPHAMERIC  TO  INTEGER  FORM 

15  NCODEKJlT( ALPHA, 1,2,1. )+0.5 

C  FIND  MAX  NUMBER  OF  DATA  ITEMS  FOR  THIS  NCODE 
DO  17  I-l.NCOm 
IF  ( ITBLE(I , 1 ) -NCODE )  17.16.17 

16  MAXNOITBLE(I.2) 

GO  TO  21 

17  CONTINUE 

C  SEARCH  DATA  ROUTINE 

C  SEE  IF  ANY  DATA  FOR  THIS  CARD 

DO  19  I-l.NCOW 
IF  (ITBLE (1,1) -NCODE)  19,18,19 

18  MAXNO-ITBLE(I.2) 

GO  TO  20 

19  CONTINUE 

20  CONTINUE 

21  IF  (MAXNO)  23.22,23 

22  RETURN 

C  ZERO  ARRAYS  AND  COUNTERS 

23  DO  47  1-1,310 

47  DATA  (I)-0. 

NDATA-1 
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2*  NCHAR-0 

25  DO  26  1-1,10 
AUXA(I)-0. 

26  AUXB(I)-0. 

IT1-1 

IT2-1 

SIGN-1. 

LDGIT-0 

KDGIT-O 

C  CARRY  OUT  DIGIT  BY  DIGIT  SEARCH  AND  ACCUMULATION 

27  NCHAR-NCHAR+1 

C  HAVE  WE  CONSIDERED  ALL  CHARACTERS  -  RETURN  IF  SO 

IF  (NCHAR-60)  26,32,1 

28  DO  29  1-1,15 

IF  (CHAR(NCBAR) -ZALPHA(I ) )  29.30,29 

29  CONTINUE 
GO  TO  32 

30  GO  TO  (33,33,33,33,33,33,33,33.33,33,32,27.36.32,31,27),  I 

C  SN  39  HANDLES  SIGN  CONTROL  ON  1130  VERSION 

31  SIGN— 1.0 
GO  TO  27 

C  CHARACTER  IS  BLANK  OR  COMA  -  DOES  IT  FOLLOW  A  DIGIT 

32  GO  TO  (27, 48).  IT1 

C  CHARACTER  IS  A  DIGIT  -  HAS  A  DECIMAL  BEEN  ENCOUNTERED 

33  GO  TO  (34,35),  IT2 

34  LDGIT-LDGIT-H 
IT1-2 

AUXA ( LDGI T ) -CHAR ( NCHAR > 

GO  TO  27 

35  KDGIT-KDGIT+1 

AUXB ( KDGI T ) -CHAR ( NCHAR ) 

GO  TO  27 

C  CHARACTER  IS  A  DECIMAL  -  DOES  IT  FOLLOW  A  DIGIT 

36  GO  TO  (37,38),  IT1 

37  IT1-2 
LDGIT-1 

38  IT2-2 

GO  TO  27 

C  ROUTINE  TO  CONVERT  ALPHABETIC  ARRAY  TO  FLOATING  POINT  NUMBER 
48  DATA  (NDATA) -GIT (AUXA , 1 , LDGIT , 1 . )+GIT ( AUXB .1.10.0.) 

DATA  ( NDATA )-DATA( NDATA )*SIGN 

C  IS  ALL  DATA  FURNISHED  YES-RETURN  NO  INCREASE  N  DATA  KEEP  ON 

IF  (NDATA-MAXNO)  41.39,39 

39  ICC-0 

40  RETURN 

41  NDATA-NDATAn 
GO  TO  25 


C 

42 

43 

44 

45 

46 


FORMAT  (2A1.9A2.60A1) 

FORMAT  ( 1H1 ) 

FORMAT  (5X.2A1.9A2.60A1) 

FORMAT  ( 2A1 , 9A2 , 60A1 ) 

FORMAT  (10X.20HCOWAND  NOT  IN  TABLE) 
END 


FUNCTION  Grr  ( TCARD . J . JLAST .SHIFT) 
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C  Convert*  alphabetic  array  to  floating  point  number 


DIMENSION  TCARD(IO),  A<10) 

DATA  ACD/1H1/  ,A(2)/1H2  /  ,A(3)/1H3  / ,  A(4)/1B4/ ,A(5>/1H5/ ,A(6) /1H6/ 
DATA  A(7)/lH7/,A(8)/lH8/fA(9)/lH9/,A(lO)/lH0/ 

GIT-0 . 

TEN- 10. 

Slft^O. 

DO  3  JNOW-J.JLAST 
TTEST-TCARD(JNCW) 

C  CHECK  FOR  LAST  ENTRY 

IF  (TTEST.EQ.O. )  GO  TO  A 
C  FIND  NUMBER  AND  COMPUTE  VALUE 

DO  2  NUMB-1 , 10 
IF  (TTEST-A(NUMB))  2,1,2 

1  Z TEST-NUMB 

IF  (ZTEST.EQ.10. )  ZTEST-0 . 

SUM-SUM* TEN+ZTEST 
GO  TO  3 

2  CONTINUE 

3  CONTINUE 

A  IF  (SHIFT)  6,5,6 

5  FI-JNOW- 1 
SUM-SUM* ( 0 . 1**FI ) 

6  GIT-SUM 
RETURN 
END 


SUBROUTINE  STHYD 

C  THIS  SUBROUTINE  STORES  THE  COORDINATES  OF  HYDROGRAPHS. 

C  ADDS  BASEFLOW  TO  HYDROGRAPH  STORED 

CCH10N/ BLOCK 1/  OCFS( 300 , 6) . DATA(31 1 ) ,CFS( 300 ) , CTBLE( 50,11), 
&RAIN(300 ) ,ROIN(6) , 

&A(20,70),Q(20,70) , DEEP (20 . 70 ) , IT8LE( 50 , 2) , DP( 20 ) , SCFS( 20 ) , C( 20 , 6 ) , 
&ZALPHA( 20 ) . I END ( 6 ),DA(6),DIST(6), SEGN(6) , DT( 6 ) , PEAK ( 6 ) , ISG(6) . 

&NPU .  NHD .  NER .  MAXNO ,  NCOW .  ICC .  NCODE  ,  TIME .  KCODE .  ICODE 

DIMENSION  DUW4Y ( 300 ) 

ID-DATA(l) 

NHD-DATAC2) 

DT(ID)-DATA(3) 

IF ( KCODE. EQ.0)GO  TO  10 
DATA(*)-DATA(A)/2. 590 
DATA( 5 )“DATA( 5 ) /0 . 02832 
DO  11  J-6,305 
DATA(J)-DATA(J)/. 02832 
11  CONTINUE 

10  DA(ID)-DATA(A) 

BSF-DATA( 5) 

C  BASEFLCW 


J-6 

C  REMAINING  DATA  ARE  FLCW  RATES 
IF(BSF.GT.O)THEN 
OCFS ( 1 , ID ) -DATA ( J ) +BSF 
GOTO  51 
ENDIF 

OCFSd.ID)-DATA(J) 

51  PEAK(ID)  -  1. 

RO  -  DATA(J) 

DO  4  1-2,300 
J-J+l 

IF(BSF.GT.0)THEN 
OCFS ( I . ID ) -DATA ( J ) +BSF 
GOTO  50 
ENDIF 

OCFS(I.ID)-DATA(J) 

50  RO  -  RO  +  OCFS (I, ID) 

C  IS  FLCW  RECEDING 

IF  (OCFS (I , ID ) -OCFS C 1*1, ID) )  1,2,2 
C  HAS  FLOW  RECEDED  TO  CUTOFF  RATE 

1  IF  (OCFS(I.ID))  5,5,4 

C  DETERMINE  PEAK  FLOW 

2  IF(0CFS( I , ID )  -  PEAK(ID))  4,4.3 

3  PEAK(ID)  -  OCFS ( I , ID ) 

4  CONTINUE 

5  IEND(ID)-I-1 
M-IEND(ID) 

ROIN(ID)  -  RO*DT ( ID ) *3600 
IF(NPU.LE . 0)GO  TO  7 
IF ( ICODE . EQ . 0 )GO  TO  6 
ROINl-ROIN(ID)*  0.02832 
DA1-DA( ID )*2 . 590 
PEAK1-PEAK( ID  )* . 02832 
DO  13  J-l.M 

DUt«Y ( J ) -OCFS ( J .  ID ) *0 . 02832 

13  CONTINUE 

WRITE( 7 . 14 ) ID , NHD . DT( ID ) . DAI . PEAK1 . ROIN1 . IEND ( ID ) , ICODE 
WRITE ( 7 , 15 ) (DUMMY (I) , 1-1  ,M) 

RETURN 

C  PUNCH  CODE 

6  WRITE(7,8)ID,NHD,DT(ID) , DA(ID) .PEAK (ID) ,ROIN( ID) . IEND ( ID) . ICODE 
WRITE  (7.9)  (OCFS(J. ID) , J-l.M) 

7  RETURN 
C 

8  FORMAT (  'RECALL  HYD‘ ,T21, 'ID-MI, T29. ’HVD  NOM3.T42, 'DT-' ,F9 
&6, '  HRS' ,T61, *DA-‘ ,F8.3, ’  SQ  MI ' /T21 , ’ PEAK- ' . F7 . 0 , 'CFS ' . T40 . ’ RO- ' 
&F6.3,"  CFS  "^59, "NO  PTS  -".I3/T21 ,  "CODE-"  .  I1/T21 . 

&"FLOW  RATES") 

9  FORMAT  (T21.7F8.0) 

14  FORMATCRECALL  HYD"  ,  T21 ,  "ID-"  .  1 1 .  T29 ,  "BYD  NO  -".I3.T42, 

&-DT-" . F9 . 6 , "HRS" , T61 , "DA-" ,F8 . 3, "SQ  KM" /T21 . "PEAK" , F7 . 2 . 
&"CMS",T40, "R0-",F6.0,"  CUMECS  ",T59."NO  PTS-" . I3/T21 . "CODE-" . 
&I1/T21, "FLOW  RATES") 

15  FORMAT  (T21.7F8.2) 

END 


SUBROUTINE  RECHD 


C  THIS  SUBROUTINE  RECALLS  PREVIOUSLY  COMPUTED  AND  PUNCHED 

C  HYDROGRAPHS 

CCfrMON / BLOCK  1  /  OCFS(300, 6)  ,DATA(311)  ,CFS(300)  ,CTBLE(50 , 11) , 
&RAIN(300),ROIN(6). 

&A( 20, 70), Q(20,70) , DEEP (20 , 70  > , ITBLE( 50 , 2 ) , DP(20 ) ,  SCFS ( 20 ) , C ( 20 , 6 ) , 
&ZALPHA( 20 ) , IEND(6) , DA( 6) ,DIST(8) , SEGN(6 ) , DT( 6 ) , PEAK( 6) , ISG( 5 ) , 

&NPU .NHD.NER, MAXNO , NCOM .  ICC , NCODE , TIME . KCODE .  ICODE 
MET1«DATA(8) 

IF (MET1 . EQ. 0 )GO  TO  2 
DATA(4)-DATA(4)/2.590 
DATA( 5) "DATA ( 5 ) / . 02832 
DATA(6)-DATA(6)/25.4 
^DATA(7> 

DO  3  I-9.M+9 
DATA( I )-DATA( I )/0 . 02832 
3  CONTINUE 

2  ID-DATA(l) 

NHD-DATA(2) 

DT(ID)-DATA(3) 

DA(ID)-DATA(4) 

PEAK(ID)-DATA(5) 

ROIN(ID)-DATA(6) 

IEND(ID)-DATA(7) 

M-IEND(ID) 

J  -  9 

C  REMAINING  DATA  ARE  FLCM  RATES 

DO  1  I-l.M 
OCFS(I,ID)-DATA(J) 

1  J-J+l 

RETURN 
END 


SUBROUTINE  CMPHYD 


C  This  subroutine  develops  a  unit  hydrograph,  converts  rainfall  data 
C  into  runoff  by  calling  the  soil  moisture  finite  difference  model, 

C  or  the  Curve  Number  routine. 

C  and  sums  these  two  to  produce  the  storm  runoff  hydrograph. 


COMON/BLOCK1/  OCFS (300 , 6 ) , DATA(31 1 ) ,CFS ( 300 ) , CTBLE( 50,11), 
&RAIN(300) ,R0IN(6) , 

&A ( 20 , 70 ) , Q(20 , 70) , DEEP ( 20 . 70) , ITBLEC  SO , 2 ) , DP(20 ) , SCFS ( 20 ) , C( 20 , 6 ) , 
&ZALPHA ( 20 ) , IEND( 6 ) , DA(6) , DIST(6 ) , SEGN(6 ) , DT(6 ) . PEAK ( 6 ) , ISG( 6 ) , 
&NPU , NHD.NER,  MAXNO . NCC**! , ICC . NCODE .TIME, KCODE , ICODE 

DIMENSION  DUT*fY<300) 

TEMP-0. 


C  Input  data  read  Into  subroutine 
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ID-DATAU) 

NHD-DATA(2) 

DT(ID)-DATA(3) 

if(KCOde.ne.o)then 

C  Convert  metric  to  imperial 

DATA(4)-DATA(4>/2.590 
IF(DATA(6) . LT. 0 )G0  TO  AO 
DATA(6)»DATA<6)/0.3048 
DATA(7)-DATA{7)/1.6 

ENDIF 

AO  DACID)-DATA(A) 

C 

CH-DATA(5) 

C  Data  items  6  and  7  normally  hold  watershed  height  and  length  and 
C  from  these  the  constants  XK(recession  constant)  and  Tp(time  to  peak) 
C  can  be  calculated. 

C  If  XK  and  Tp  are  known  however,  they  can  be  entered  instead 
C  and  a  negative  sign  is  put  before  their  values. 

IF  (DATA(6) . LT.O. )THEN 
XK— DATAC6) 

TP—  DATA(7) 

ELSE 

HT-DATA<6) 

XL-DATA(7) 

SLOPE-HT/XL 

XLDW-(XL**2.)/DA(ID) 

XK-27 .0*(DA(ID)**. 23  D* (SLOPE** (- . 777 ) )* (XLDW** . 12A ) 

TP-4 . 63*(DA( ID )** . *22)* (SLOPE** ( - . A6 ) )*(XLDW** . 133 ) 

ENDIF 


C  The  storm  runoff  array  is  intialised  to  0,  and  peak  of  hydrograph  to  1 

DO  A  1-1,300 
A  OCFSCI , ID)“Q . 

PEAK(ID)«1. 

C  Compute  '  N*  by  iteration 
XN-5.0 
XKTP-XK/TP 
DO  6  1-1,50 

TINF-1 . +SQRT< 1 . / (XN-1 . ) ) 

XN1-. 05/<XKTP*(ALOG(7ZNF/(TINF+. 05) )♦.  05) )  +  2 . 

DIFF-ABS(XNl-XN) 

IF  (DIFF-.001)  7,7.5 

5  XN-XN1 

6  CONTINUE 
WRITE  (6.29) 

29  FORMAT ( *  N  DID  NOT  CONVERGE  AFTER  50  ITERATIONS.’) 

GO  TO  28 
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C  Compute  ’Cl' 

7  DELT-TINF/100. 

TC1-0 . 

XN1P-XH-1. 

XN1W-1.-XN 

DO  a  1-2,101 
TC1-TC1+DELT 

8  CFS(I>-(TC1**XN1P)*EXP(XN1M*(TC1-1.  )) 

SUM-CTS(101)/2. 

DO  9  1-2.100 

9  SU^SUM+CFSU) 

C1-SUM*DELT 

C 

C  Compute  *  B ’ 

CFSII-CFSdOl) 

TTINF-TINF*TP 
TREC1-TTINF+2 . *XK 
EEE-EXP<  (TTINt  ^C1)/XK) 

XK1-3 . *XK 

B-645 . 333/ (C1+CFSII*(XKTP* ( 1 . -EEE)+EEE* (XK1/TP) ) ) 

C 

C  Compute  'QP’  and  ‘CFSI’ 

C 

QP-(B*DA(ID))/TP 
CFSI-QP*CFS ( 101 ) 

CFSR1-CFSI*EEE 
IF ( ICODE . EQ . 0 )GO  TO  45 
QPl^JP* . 02832 
WRITE(6 , 38)XN,QP1 

38  FORMAT  ( '  Shape  constant,  N  -  ’^6.3/’  Unit  peak  -  ’.F10.1.1X 

& , ' cms ' / ) 

GO  TO  44 

45  WRITE  (6,30)  XN.QP 

30  FORMATC  Shape  constant,  N  -  ’  ,  F6 . 3 /  *  Unit  peak  -  '.F10.1.1X 
* , ' cms ’ / ) 

C 

44  CONTINUE 

C 

C  Determine  the  incremental  runoff 
C 

IF (KCODE . NE . 0 )THEN 

IF(DATA(8) . LT. 0)GO  TO  13 

C  Convert  rainfall  data  from  mm  to  inches. 

DO  34  K-8,308 

DATA(K)«DATA(K)/25.4 

34  CONTINUE 
ENDIF 

C 

35  J-8 

IF  (OATA(J))  13,10,10 

10  RAINd)-DATA(J) 

DO  11  1-2,300 

J-J+l 

RAIN(I)-DATA(J) 

IF  (RAIN(I)-RAIN(I- 1) )  12.11.11 

11  CONTINUE 

12 
C 


NUMB- I - 1 
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l 


i 

I 

\ 


I 


i 

r 


c 

C  Curve  number  routine 

13  IF (CN . LE . 0  >GOTO  2Q1 

C  STORAGE 

R-1000./CN-10 
B1-.2*R 
DO  15  1*1 , NUMB 
IF(RAIN( I )"B1 )33 ,33,14 
33  DATA(I)-0 

01-0 
GOTO  15 

14  02-((RAIN(I)-Bl)**2.  )/(RAIK(I)-f.8*R) 

DATA(I)^32-Q1 

Ql-02 

15  CONTINUE 
GOTO  202 

C 

c 

201  DO  5555  1-1,300 

5555  DATA(l)-0 

TEMP-DT(ID) 

C 

CALL  SOILM( TEMP, NUMB. RAIN, DATA) 

C  If  no  runoff  has  been  generated  then  the  simulation 
C  stops. 


DO  100  1*1 , NUMB 
IF(DATA( I ) . EQ. 0 . )GOTO  100 
GOTO  200 
100  CONTINUE 

WRITE t 6, 300) 

300  FORMAT (’  model  generated  no  runoff'/ 
&'  Simulation  terminates') 

STOP 

200  CONTINUE 
Compute  unit  hydrograph 


202  T2-0. 

CFS( 1 )*0 . 

DO  20  1-2,300 
T2-T2+DT(ID) 

IF  (T2-TTINF)  16, 16, 17 

16  CFS(I)-QP*((T2/TP)**XN1P)*EXP(XN1M*(T2/TP-1. )) 

GO  TO  20 

17  IF  (T2-TREC1)  18,10,19 

18  CFS( I )“CFSI *EXP( (TTINF-T2) /XK) 

GO  TO  20 

19  CFS(I)-CFSR1*EXP((TREC1-T2)/XK1) 

IF  (CFS(I)-l.)  21,21,20 

20  CONTINUE 
1-300 

21  ICND-I 
C 

C 

C  Compute  the  storm  runoff  hydrograph  by  summing  the  unit  hydrograph  and 
C  the  runoff. 


a 
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c 

c 

DO  24  J-2 , NUMB 

N-J+ICHD-2 

IF  (N-300)  23.23.22 

22  N-300 

23  1-2 

DO  2*  K-  J.H 

OCFS (K , ID )  -OCFS (K , ID ) +DATA ( J ) *CFS ( I ) 

I-I+l 

24  CONTINUE 
C 

C  Compute  the  runoff  volume  and  determine  the  peak. 

C 

C 

RO  -  0. 

DO  26  I  -  2 , N 

RO  -  RO  +  OCFS (I , ID) 

IF  (OCFS( I , ID)-PEAK( ID ) )26 , 26 , 25 

25  PEAK(ID)-  OCFS ( I , ID) 

26  CONTINUE 

I END  (ID)  -  N 
ROIN(ID)«RO*DT(ID)*3600 
C 

C  PUNCH  CODE 

IF  (NPU)  28,28.27 

27  IF ( ICODE . EQ . 0 )GO  TO  39 
ROIN1-ROIN ( ID  >*0 . 02832 
DA1-DA( ID >*2.590 
PEAK1-PEAK (ID )*. 02832 

DO  41  J-l.N 

DU^(  J )-OCFS(  I .  ID)*0 . 02832 
41  CONTINUE 

WRITE ( 7 , 37 ) ID . NHD , DT( ID ) . DAI . PEAK1 , R0IN1 , IEND( ID ) , ICODE 
WRITE (7 ,42) (DUM1Y( I ) . 1-1 , N) 

RETURN 

39  WRITE(7,31)ID,NHD,DT(ID),DA(ID).PEAK(ID).R0IN(ID),IEND(ID).IC0DE 

WRITE  (7,32)  (OCFS(I.ID) ,I-1,N) 

28  RETURN 


C 

31 


37 


42 

32 


FORMAT (  ’RECALL  HYD’  ,T21,  '  ID-'  ,  I1.T29.  ’HYD  NO’  ,  13 ,  T42 ,  ’  DT-’  .  F9. 
&6,'  HRS’ ,T61, ’DA-' ,F8.3, '  SQ  MI • /T21 . *  PEAK-' , F7 . 0 . 'CFS ’ . T40 . ’ RO- ’ . 
&F6.3,'  CFS'.TSO.'NO  PTS-'.I3/T21 ,  'CODE-" , I1/T21 FLOW  RATES') 
FORMAT (  ’RECALL  HYD ' . T21 . ' ID-’ . II , T29, ' HYD  NO ’ , 13 . TA2 . ’ DT- ’ . F9 . 
&6.’  HRS ’ , T6 1 , ’ DA- ’ , F8 . 3 . '  SQ  KM’ /T21. ’ PEAK-’ ,F7 . 2, ’CMS  ’ , T*0 , -RO-’ . 
&F6.0.’  CUMECS  \T59.’NO  PTS*'  ,  I3/T21.  *CODE--’  .I1/T21 .  'FLOW  RATES’) 
FORMAT  (T21.7F8.2) 

FORMAT  (T21.7F8.0) 

END 


SUBROUTINE  SOILM(DT, IR .CUMRAIN, DATA) 

C  A  phy.lc.lly  b.,.d  p.r«m.t.r  .ntiUr.tion  mod.:  which  .imuUt.,  n..r  ,urf.c 
c  *oil  water  movement,  and  hence  runoff. 


* 


C  Variables  used  in  this  subroutine 


C  TIME 

C  SRI 

C  SR2 

C  SR3 

C  NLA 

C  MLB 

C  ML 

C  SATCON 

C  SATCON2 

C  SATCON 3 

C  EMAX 

C  SIMDUR 

C  DETCAP 

C  AF 

C  WT 

C  THETA 

C  TCOM 

C  ALR 

C  AMR 

C  NQ 

c  x 

C  Y 

C  X2 

C  Y2 

C  X3 

C  Y3 

C  IR 

C  DT 

C  CUMRAIN 

C  NSCOL 

C  I PC AREA 

C  IOUT 

C 
C 


Time  whan  simulation  begins  (hours). 

Soil  water  content  at  saturation  layer  1. 

(m3/m3)  layer  2. 

layer  3. 

Number  of  cells  in  layer  1. 

Number  of  cells  in  layer  2. 

Total  number  of  cells  in  coluosi 
Saturated  permeability  (ms-l)  layer  1. 

layer  2. 
layer  3. 

Maximum  evaporation  during  the  day  (ms-1). 

Simulation  duration  (hours). 

Surface  detention  capacity  (m) . 

Simulation  iteration  period  (secs). 

Write-out  time  period  (hrs). 

Initial  soil  water  content  for  each  cell  (ro3/m3). 
Thickness  of  each  cell. 

Rain  start  time  (hours). 

Rain  stop  time. 

Number  of  observations  on  suction  moisture  curve. 

Moisture  values ....  layer  1  (m3/m3). 

Suction  values. ....  layer  1  (bars). 

layer  2. 
layer  2. 
layer  3. 
layer  3. 

Number  of  rainfall  observations. 

Rainfall  data  time  increments  (hours). 

Cumulative  rainfall  data  at  DT  time  increments  (inches). 
Number  of  soil  columns. 

Percent  area  of  soil  column. 

Determines  amount  of  output. 

1  -  total  output 
0  -  shorter 


C  Note: 

C  If  SRI,  SR2.  SR3,  SATCON,  SATC0N2 ,  SATCON 3 ,  DETCAP.  THETA,  X,  X2.  Or  X3 
C  are  proceeded  by  an  ’A' ,  then  the  variable  type  is  double  precision 
C  rather  than  real.  If  SRI,  SR2,  SR3.  SATCON.  SATC0N2 ,  SATC0N2 ,  DETCAP, 

C  OR  THETA  are  preceeded  by  an  'S’, then  the  variable  represents  the 

C  standard  deviation  of  that  particular  soil  hydrological  characteristic. 


C  SCURV1  Standard  deviation  of  soil  moisture  curve  for  layer  1 
C  SCURV2  laYer  2 
C  SCURV3  IaY®r  3 


INITIAL  SECTION 
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c 

DIMENSION  FLUX (20) , TCCM(20 ) , SWP(20 ) , THETA (20  > ,COND(20 ) 
DIMENSION  VOL ( 20 ) , ANFLUX ( 20 ) , AVCOND ( 20 ) , DEPTH <  20 ) , DIST ( 20 ) 
DIMENSION  X(20),Y(20),G(20),GZ(20), FSWP(20 ) , CNT(20 ) 
DIMENSION  CUMRAIN(251),Z(20),PPT(250)IXP(20).FS(20) 
DIMENSION  DATA(300),WDATA(300,10),BPOT(20) 

DIMENSION  G2(20 ) , Y2(20) ,X2(20 ) , G22(20 ), 22(20 ) 

DIMENSION  G3 ( 20 ) , Y3 ( 20 ) , X3 ( 20 ) , GZ3 ( 20 ) , Z3 ( 20 ) 

DIMENSION  RSAT ( 20 ) 

DIMENSION  AX (20) ,AX2(20) ,AX3(20) ,ATHETA(20) 

DIMENSION  XNEW( 20 ) , YNEW(20 ) , X2NEW( 20 ) , Y2NEW( 20 ) , 

&  X3NEW(20) ,Y3NEW(20) 

DOUBLE  PRECISION  G05DDF 
DOUBLE  PRECISION  D LOG 10 

DOUBLE  PRECISION  ATHETA , AX , AX2 , AX3 , ADETCAP , ASR1 , ASR2 , ASR3 , 

*  ASATCON , ASATCON2 , ASATCON3 , BSATCON , BSATCON2 , BSATCON3 . 

*  SDETCAP , SSR1 , SSR2 . SSR3 , STHETA , SSATCON , SSATCON2 , SSATCON3 , 

*  SCURVl , SCURV2 , SCURV3 
C 

C 

C 

C  READ  IN  DATA 

C  . 

C 

C 

C 

READ (2 5, * 'TIME , ALR, AMR, SIMDUR 
R£AD(25,*)IOUT 
READ(25,*)AF.WT 
READ(25,*)NSCOL 


C  The  array  RAIN  which  is  passed  to  the  subroutine  as  a  cumulative 
C  rainfall  total  is  in  inches.  This  has  to  be  transferee!  to  array 
C  PPT  which  is  in  m  and  represents  the  total  for  each  time  increment. 
IRR-IR-1 
DO  100  I-l.IRR 

100  PPT(I)»(CUMRAlN(in)-CUMRAlN(I))V025* 

DO  3*543  W-l.NSCOL 

C  For  each  soil  column  in  turn,  read  in  data  and  proceed  through 

C  simulation  to  determine  runoff 


READ (25 , * ) I PC AREA 
READ(25,*)NL,NLA,NLB 
READ (25, *) (TCOM(I) , I*1,NL) 

READ ( 25 . * ) EMAX . ADETCAP , SDETCAP 

READ (25 , ♦ ) ASR1 . SSR1 , ASR2 , SSR2 , ASR3 , SSR3 

READ ( 25 , * ) ASATCON , SSATCON . ASATCON2 , SSATCON2 , ASATCON3 , S^ATCON3 
READ(25 , * > (ATHETA( I ) , 1-1 . NL ) 

READ (2 5,*) STHETA 
R£AD(25, *)NQ 
R£AD(25.*)(AX(I),I«1.NQ) 

READ (25 ,#)(Y(I),I*1,NQ) 

READ (25 , * ) SCURVl 


READ (25 . *) (AX2 ( I ) , 1-1 , NQ) 
R£AD(25,*)(Y2CI).I-1,NQ) 
READ ( 25 , * )SCURV2 
READC25.*)(AX3(I).I-l(NQ) 
READ (25, *)(Y3(I) ,I«1 (NQ) 
READ(25.*)SCURV3 

NQJ-NQ 

NLL-NL+1 

IF(AMR.LT.ALR)THEN 
AMR-AMR+24 . 0 

END  IF 

C 

C 

c  CHECK  DATA  INPUTS 

C  _ 

C 

NERROR'O 


u  ‘-beck  number  of  cells  in  soil  column 
IF ( NLA+NLB . GE . NL ) THEN 
WRITE(6, 1015) 

1015  FORMAT ( '  Error-NLA , NLB ,  NL ' ) 

NERROR-NERROR+1 


C  Check  dimensions  of  input  vectors 

IF<NO.GT.20 .OR. NL .GT.20 .OR. IR.GT. 250) THEN 
WRITE(6, 1020) 

1020  FORMAT ( *  Error-limit  exceeded , NQ.NL,  IR'  ) 
NERROR-NERROR+1 


C  Check  rainfell  passed  from  CMPHYD 
kn-ir-i 
DO  50  1*1,  KN 

IF(CUMRaIN( 1+1) . LT.CUMRAIN(I) )THEN 
WRITE(6, 1030) 

1030  FOSHATC  Error-not  c«ul,uv,  rainfall  toils') 

NERSOR-NERRORn 
ENDIF 

30  CONTINUE 


C  CWk  that  initial  moiatur.  cont.nt  of  a.ch  call  li.a  -ithin  the  range  of 
C  th.  auction  moiatur.  curv,  and  do.a  not  .read  at.t.d  a.tur.t.d  Mature 
C  content. 

DO  51  1*1 , NLA 

IF(ATHETAil)  .GT.ASRDiHEN 
WRITE<6. 1050) 

1050  FORMAT ( ’  Error-THETA  larger  then  sat  moisture  content(l)-) 

NERROR-NERROR+1 
ENDIF 

IF  (ATHETA(I) .GT.AX(NQ) . OR  . ATffETAC I ) . LT . AX( 1 ) ) THEN 
WRITEC6. 1055) 

1055  FORMAT ( '  Error-THETA  outside  range  of  curves-d)’) 

ENDIF 
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T  " 


51  CONTINUE 
NLAA-NLA+1 
NLH-NLA+NLB 

DO  52  I-NLAA , NLH 

IF < ATHETA ( I ).GT.ASR2) THEN 
WRITE (6, 1060) 

1060  FORMAT ('  Error-THETA  larger  than  sat  moisture  contant(2)') 

NERRGR-  NEKROR + 1 
END  IF 

IF ( ATHETA (I ). GT . AX2 ( NQ ). OR . ATHETA ( I ). LT . AX2 ( 1 )) THEN 
WRITE (6, 1065) 

1065  FORMAT ('  Error-THETA  outsida  rang*  of  curve- (2)') 
nerror-nerror+i 
ENDIF 

52  CONTINUE 
NLBB-NLB+NLA+1 
DO  53  I-NLBB.NL 

IF ( ATHETA ( I ). GT . ASR3 ) THEN 
WRITE (6, 1070) 

1070  FORMAT ( '  Error-THETA  larger  than  sat  moisture  content(3)’) 

STOP 
ENDIF 

IF ( ATHETA ( I ) . GT . AX3 ( NQ ) . OR . ATHETA (I ) . LT . AX3 ( 1 ) ) THEN 
WRITE(6, 1075) 

1075  FORMAT('  Error-THETA  outside  range  of  curve  -(2)’) 
NERROR-NERROR+1 

ENDIF 

53  CONTINUE 
C 

IF  (NERROR.NE.O)THEN 
WRITE(6,1076)NERROR 

1076  FORMAT ( ’  SOILM:  number  of  input  data  errors  *,I2, 

&' Simulation  terminates') 

STOP 

ENDIF 

C 

C 

C 

C 

C  DEPTH  CALCULATION 

C  - . . 

C 

C 

C 

C  The  variable  DEPTH  is  calculated.  This  refers  to  the  distance  from 
C  ground  level  to  any  cell  midpoint. 

C  DIST  refers  to  the  distance  between  any  two  adjacent  cell  midpoints. 
C 

DIST(l)«TC0M(l)/2. 

DEPTH (l)-DIST(l) 

DO  110  I-2.NL 

DEPTH ( I ) -DEPTH (I - 1 >+0 . 5* ( TC0M< I - 1 ) ♦TCCM( I ) ) 

110  DIST< I )-0 . 5*(TC0M(I-1)+TC0M(I ) ) 

C 

C 

C 

C 

C  PARAMETER  VARIABILITY 


JOk _ _ 


C  - - 

c 

C  Five  input  variables,  detention  capacity,  soil  water  content  at 
C  saturation,  soil  moisture  content  at  given  tensions,  saturated  conductivity 
C  and  initial  moisture  content  are  varied  stochastically. 

C  NAG  functions  are  called  which  return  a  ’psuedo  random'  value  from  a 
C  distribution  with  a  given  standard  deviation  and  mean. 

C  All  are  assumed  to  have  a  normal  distribution  except  the  seturated 
C  conductivity  which  takes  on  a  lognormal. 

C 

C 

c 

C  Generate  only  one  set  of  stochastic  variahlas  to  run  in  HYMO. 

C 

C 

C  RANDOM  PARAMETER  VALUE 

C  . . - . 

C 

WRITE (6, 1079) 

1079  FORMAT < *  INCREMENTAL  RUNOFF-Parameter  variability  included’ // ) 

C 

C  Detention  capacity. 

DETCAP-GO  5DDF ( ADETCAP , SDETCAP ) 

IF ( DETCAP . LT . 0 . ) DETCAP- 0 . 0 

SD-SDETCAP 

WRITE (6,1 180 )SD 

1180  FORMAT ('  SD  of  detcap  \F5.3) 

C 

C  Soil  water  content  at  saturation 

SRI-G05DDF(ASR1,SSR1) 

SR2-G0 5DDF ( ASR2 . SSR2 ) 

SR3-G0  5DDF ( ASR3 , SSR3 ) 

SD1-SSR1 

S02-SSR2 

SD3-SSR3 

WRITE(6, 1181)SDl,SD2,SD3 

1181  FQRMAT('  SD  of  saturated  soil  content ’ ,F5 . 3 , ’  layer  1’/ 

&  ’  *  , F5.3, ’  layer  2’/ 

&  '  ' ,F5.3 , ’  layer  3’ ) 

C 

C  Soil  moisture  content  at  given  tensions 

C  Layer  1 

CALL  SMCURV(SR1 ,NQ, AX, Y.XNEW, YNEW.SCURV1) 

DO  120  1-1,20 
XU)-XNEW(I) 

120  Y(I)-YNEW(I) 

C  Layer  2 

CALL  SMCURV ( SR2 ,NQ, AX2.Y2 . X2NEW , Y2NEW , SCURV2 ) 

DO  130  1-1,20 
X2( I )*X2NEW( I ) 

130  Y2( I )-Y2NEW( I ) 

C  Layer  3 

CALL  SMCURV ( SR3 , NQ , AX3 . Y3 , X3NEW . Y3NEW , SCURV3 ) 

DO  HO  1-1.20 
X3 ( I )-X3NEW( I ) 

140  Y3(I)-Y3NEW(I) 

SD1-SCURV1 

SD2-SCURV2 


SD3-SCURV3 

WRITE(6, 1182)SDl,SD2,SD3 

1182  FORMAT ('  SD  of  suction  moisture  curve’,  F5.3,'  layer  1'/ 
&  *  ' ,  F5.3.'  layer  2*/ 

&  '  ’ ,F5.3, '  layer  3' ) 

C 

C  Saturated  conductivity  for  each  layer 
BSATCON-DLOGIO ( ASATCOH ) 

SATCON-G05DDF ( BSATCON , SSATCON ) 

SATCON-10**SATCON 
BSATCOH  2»DLOG 1 0 ( ASATCON 2 ) 

SaTCON2 -CO 5DDF (BSATCOH* , SSVTCON2) 

SATCON2- 1 0 * *SATCON2 
BSATCON3-DLOG10 ( ASATCON3 ) 

SATCON3-GO  5DDF ( BSATCOH 3 , SSATCON 3 ) 

SATCON3-10**SATCOH3 
SD1-SSATC0H 
SD2*SSATCON2 
SD  3  "•SSATCON  3 

WRITERS, 1183 )SD1 , SD2 , SD3 

1103  FORMAT ( ’  SD  of  sat  conductivity* ,F5. 3. •  lay^r  V  t 

&  ’  ’ ,F5. 3, ’  layer  2'/ 

&  *  ' ,F5.3, ’  layer  3* ) 

C 

C  Initial  moisture  content 
DO  150  I-l.NL 

150  THETA  ( I )  «<50  5DDF  ( ATHETA  ( I ) .  STHET A ) 

C  Check  on  initial  soil  moisture  values 

DO  160  I-l.NLA 

IF(THETA(I).GE.X(20))THETA(I)-XC20)-0.001 
160  IF ( THETA ( I ) . LE.XC 1 ) )TRETA( I )*Xf 1 j+0 . 001 

DO  170  I-NLAA.NLH 

IF(THETA(I).GE.X2(20))TBETaCI)*X2(20) -0.001 
170  IF ( THETA( I ) . LE . X2 ( 1 ) ) THETA( I >“X2 ( 1 >  +0 . 00 1 

DO  100  I-NLBB.NL 

IF(THETA(I) .GE.X3(20) )THETa(I)“X3(20)-0. 001 
180  IF(THETA(n.LE.X3(l))THETA(I)-X3(l)+0.001 

SD-STHETA 
WRITER,  11S*)SD 

1184  FORMAT ( 1  SD  of  initial  water  content 1 . F5 . 3 ) 

C 

C 

C 

C  HYDRAULIC  CONDUCTIVITY  CALCULATION 

C  . . . 

c 

c 

c 

C  The  hydraulic  conductivity  is  calculated  from  surkion  moisture 
C  data  tor  each  layer. 

NQJ-NQ 

CALL  HYDCON(X , SATCON , SRI , Z , Y ) 

CALL  HYDCON ( X2 , SATCON2 , SR2 . 22 , Y2 ) 

CALL  HYDCON (X3.SATCON3.SR3. 23. Y3) 

C 

C 

C 

c 
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C  WRITE-OUT  INITIAL  CONDITIONS 

C  . . - . 

c 

c 

c 

C  Write-out  suction  moisture  curve  and  generated  K-vaiuea. 

C 

WRITE/6, 1080) 

1080  FORMAT ( 1 OGENERATED  K-MOISTURE  CURVE ’ / 

&'  Mil ling ton -Quirk  Method'/ 

&'  Layer  1’ , 26X, 'Layer  2' , 26X , ' Layer  3*/ 

&3('  Moisture  Suction  Unsat  K  ')) 

DO  175  1-1.20 

175  WRITE(6.1090)X(I).Y(I),Z(I).X2(I),Y2(I),22CI).X3(X).Y3(I).Z3(I) 

1090  FORMAT ( 1H  , 3(F6 . 3 . 2X , F8 . 3 , F15. 12, 2X) ) 

C  Write-out  start  conditions. 

C 

WRITE (6, 1100) 

1100  FORMAT ( 'OSTART  CONDITIONS  '/) 

WRITE(6, 1110 )TIME 

1110  FORMAT ( *  Simulation  start  time' . FA . 1 , 'hrs ’ ) 

WRITE(6. 1130 )ALR,  AMR 

1130  FORMAT/*  Precipitation  begins  at  ’ ,FA. 1.2X. ’and  ends  at  ’,FA.l) 
WRITE(6, 1140 )DT 

11A0  FORMAT ( '  Rainfall  data  time  increment  -  * ,F6. A ,2X, 'hrs ’ ) 

WRITE(6, 1120)AF 

1120  FORMAT (*  Time  increment  for  iteration  period  -  ’.F6.1, 

&2X,  ’ secs'/) 

WRITE  C  6 , 1 1 50 ) EMAX . DETCAP 

1150  FORMAT ( ’  Maximum  evaporation  during  the  day  -  ’  F10 . 8 , 2X . ’ ms- 1 ' / 

&’  Surface  detention  capacity  -  ' 1F6.A12X,’m’//> 

C 

C  Calculate  Initial  relative  saturation  of  each  cell  in  soil  column 
DO  1151  I-l.NL 

IF(X.LE.NLA)RSAT(I)-THETA<I>/SR1 

IF( I . GT . NLA . AND . I . LT . NLBB )RSAT( I )-THETA( I ) /SR2 

I F (I . GE . NLBB )RSAT ( I ) -THETA ( I ) / SR3 

1151  CONTINUE 

WRITE(6. 1152) 

1152  FORMAT ('  INITIAL  SOIL  COLUMN  CONDITIONS*//) 

WRITE ( 6 , 1153 ) 

1153  FORMAT < i IX , ’ SAT ' , 8X , ’ SAT  HYD’ , 6X. 'CELL* , IX, 'DEPTH' . 

&2X, ' INITAL' . 2X, 'REL’/ 

&1H  , 10X, ’THETA* ,7X, 'COND’ ,9X, 'NO' . 10X, ’THETA’ ,2X, 'SAT' / 

&1H  . 10X, 'm3/i»3' . 7X. 'ms-1 ’ . 1AX, ’ ro* , 5X, '»3/m3 ’ / ) 

WRITE(6. 115A)SRl.SATCON,DEPTH(l) .THETA(l) ,RSAT(1) 

U5A  FORMAT  ( ’  Layer  1  ’  ,  F7  .  A  ,  IX ,  F15. 12, 3X,  *  1 '  ,  2X ,  F6  .  A  ,  IX  .  F7  .  ft  ,  IX  ,  F5  3) 

IF(NLA.GT.  DTHEN 
DO  1155  1-2 , NLA 

WRITE ( 6 , 1156)1, DEPTH ( I ) , THETA ( I ) , RSAT ( I ) 

1156  FORMAT ( 1H  . 3AX . 12 . 2X , F6 . A . IX . F7 . 4 . IX . F5 . 3 ) 

1155  CONTINUE 

END  IF 

WRI TE  ( 6 , 1 1 57 >SR2 . SATCCN2 . NLAA .DEPTH ( NLAA ) , THETA ( NLAA ) . RSAT { NLAA ) 

1157  FORMAT ( '  Layer  2  ' ,F7. 4 , 1X.F15. 12 . 2X , 12. 2X.F6 . 4 , IX. F7 1X.F5.3) 

IF (NLB.GT.  DTHEN 

DO  1158  I-NLA+2 . NLH 


WRITE(6, 1159)1 ,DEPTH(I) .THETA (I ) ,RSAT(I) 

1X59  FORMAT ( 1H  . 34X, 12, 2X.F6. 4 , IX, F7 . 4 , IX, F5. 3 ) 

1158  CONTINUE 

ENDIF 

WRITE ( 6 . 1 1 60 ) SR3 , SATCON3 , NLH+1 . DEPTH ( NLH+ 1 ) , THETA ( NLB+ 1 ) , 

&RSAT( NLH+1) 

1160  FORMAT C  Layer  3  ' ,F7. 4 . 1X.F15. 12, 2X, 12. 2X.F6. 4 . IX, F7 . 4 , IX, F5. 3) 

IF ( ( NL-NLH ) . GT . 1 ) THEN 

DO  1161  I-NLH+2.NL 

WRITE(6, 1162)1 , DEPTH (I ) ,THETA(I) ,RSAT(I) 

1162  FORMAT ( 1H  , 34X , 12 , 2X , F6 . 4 . 1X.F7 . 4 , IX , F5. 3 } 

1161  CONTINUE 
ENDIF 

C 

C 

C 

C  INITIALISATION  OF  VARIABLES 

C  . . 

C 

C 

C 

DO  184  >1,300 

184  iwww-W 
WDATA(I,iwww)-0.0 

WATI-0 . 0 
h*M-2 

DO  185  I-2.NL 

185  ANFLUX ( I )-0 . 0 
CTIME-TIME*3600 
SRAIN 1-0.0 
CUMDRN-0 . 

CINFIL-0. 

SUMD-0. 

ICOUNT  -0 
BR-AMR-ALR 
EVAPI-0 . 0 
SOG-THETA  Cl)/ SR 1 
RTOT-O . 0 
ANFILT-0 . 0 
rrTT-0 . 0 
TG-0.0 

C 

C 

C 

C  BALANCE  CHECK 

C  . 

c 

c 

C  A  calculation  for  the  water  balance  check. 

C  The  initial  soil  water  content  of  the  soil  column. 

C 

DO  190  I-l.NL 

190  WATI-TCOM(I)*THETA<I)*WATI 
C 

c 

r 

C  CURVE  GRADIENTS 
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C  Calculations  of  the  gradients  of  the  suction-moisture  curve  and  the 
C  K-moisture  curve  for  each  layer. 

C 

CALL  GRAD(G,GZ,Y,X.Z) 

CALL  GRAD (G2 , GZ2 , Y2 , X2 , Z2 ) 

CALL  GRAD(G3 , GZ3 , Y3 ,X3 , Z3 ) 


DYNAMIC  SECTION  -  SIMULATION 


C  This  loop  is  completed  for  each  time  increment  until  end  of  simulation. 
C 

ITMAX-SIMDUR*3600/AF 
DO  9995  II-l.ITMAX 
ICOUNT-ICOUNT+AF 


CALCULATE  WATER  VOLUME  OF  EACH  CELL 


DO  200  I-l.NL 

VOL  ( I )  -TCOM ( I )  *THETA  ( I ) 


2 A -HOUR  CLOCK 


C  Calculate  REAL  TIME  for  current  iteration  period  using  the  24-hour  clock 
C 

CTIME“CTIME+AF 
IF  (CTIME.GE. 86400 )THEN 
CTIME-CTIME-86400 


SWP.HPOT.COND  CALCULATIONS 


C  Calculate  the  soil  water  pressure,  hydraulic  potential  and  conductivity 
C  for  each  cell  as  conditions  change  during  the  simulation. 


CALL  TWO(l. NLA. THETA. X.SWP.Y.G.HPOT. DEPTH. GZ.COND.Z) 
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CALL  TWO (NLAA.NLB, THETA ,X2 , SWP , Y2 ,G2 ,HPOT , DEPTH , GZ2 . COND . Z2 > 

CALL  TWO(NLBB.NL, THETA. X3, SWP, Y3,G3, SPOT, DEPTH, GZ3, COND. Z3) 

C 

C 

C 

c 

C  DETERMINE  RAINFALL 

C  - 

c 

c 

c 

C  Determine  rainfall  per  second  at  end  of  the  current  iteration 
C  period. 

C  T1  is  the  time  in  hours  when  the  current  iteration  period  ends. 

C  Check  that  T1  is  between  the  rain  start  and  atop. 

C  If  it  is,  decide  which  element  of  PPT  array  the  data  is  to  be  taken  from 
C  and  make  SRAIN  equal  to  that  precipitation  per  second. 

C  If  it  is  not  within  the  storm  period,  set  SRAIN  to  0. 

C 

C 

T1-T*AF/3600.0 

IFCTl.LE. (ALR-TIME) .0R.T1.GT. ( AMR - T IME )> THEN 
SRAIN-O . 0 
ELSE 

T2-TWAF/3600.) 

IELEM- (<T2-( ALR-TIME ))/DT)-n 
SRAIN-PPTi IELEM) / <DT*3600 . 0 ) 

END  IF 
c 
c 

C  Increment  precipitation  total  by  amount  of  precipitation  in  current 
C  iteration  period. 

C 

PPTT-PPTT+ ( SRAIN* AF ) 

C 

c 

c 

C  AVERAGE  HYDRAULIC  CONDUCTIVITY 

C  - - - - - 

c 

c 

c 

C  Average  hyraulic  conductivity  for  flow  through  boundary  between 
C  adjoining  calls  is  weighted  according  to  its  thickness. 

C 

DO  210  I-2.NL 

210  AVCOND ( I )- (COND  (1-1 ) -TC0MC  I- 1  }-KJOND(  I  )  *TCCM(  I  > ) 

&/ <  TCOM ( I - 1 ) +TC0M(  I ) ) 

C 

c 

c 

C  BOTTOM  BOUNDARY  COND IT ION 

C  . - . 

c 

c 

C  Determine  the  bottom  boundary  condition  under  the  assumption  that 
C  water  is  flowing  out  of  the  soil  column  under  gravity. 

C 


▼ 
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FLUX ( NLL )  -COND ( NL ) 

C 

c 

c 

C  FLUX  BETWEEN  CELLS 

C  - 

c 

c 

c 

C  The  flux  between  aach  call  than  follows  Darcy's  law  in  dlscrata  form. 
C 

DO  220  I-2.NL 

220  FLUX(I)-(HPOT(I-l)-HPOT(I) )*AVCOND(I)/DIST(I ) 

C 

c 

c 

C  DETERMINE  TOP  BOUNDARY  CONDITIONS 

C  - - — . 

c 

c 

c 

C  Calculate  tha  infiltration  capacity. 

C 

BNCAP- (0.0-HPOT(1))*0.5*<  S ATCON+COND < 1 ) ) / D I ST { 1 ) 

C 

C  Calculate  precipitation  excess 
C 

IF ( SRAIN1 . EQ . SRAIN ) THEN 

SUMD- < SRAIN- ANFILT ) *AF+SUMD 
ELSE 

SUMD-0 . 0+SUMD 
ENDIF 

SRAIN 1 -SRAIN 
C 

C  Calculate  amount  detained  on  the  surface. 

C 

IF (SUMD .  LT . 0 . 0)THEN 
DET AIN-0.0 
ELSE 

DETAIN-SUMD 

ENDIF 

C 

C  Calculate  evaporation,  the  flux  into  cell  1  and  runoff. 

C 

IF (SRAIN . GT . 0 . 0 )  THEN 
C 

EVAP-  0.0 
C 

IF ( SRAIN . LT . BNCAP . AND . DETAIN . LE . 0 . 0 )THEN 
ANFILT-SRAIN 
ELSE 

ANFILT-BNCAP 

ENDIF 

FLUX( 1 J-ANFILT 
C 

IF ( DETAIN .GT  DETCAP)THEN 
SUMD-DETCAP 
DETAIN-DETCAP 


/ _ a 


RUNOFF-O . 0 

IF (SRAIN . GT . BNCAP  J RUNOFF- ( SRAIN-BNCAP ) *AF 
RTOT-RTOT+RUNOFF 
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ELSE 

RUNOFF* 0 . 0 
ENDIF 
C 

ELSE 

C 

RUNOFF-O . 0 
C 

IF(CTIME.GT. 64800. AND. CTIME.LE. 21600) THEN 
EVAP-EMAX/100. 

ELSE 

EVAP-EWAX*SIN(2.*3.14159*(CTIME-21600. >/86400. ) 

ENDIF 

C 

IF (DETAIN. IE. 0.) THEN 
ANFILT-0.0 
FLUX(1)-EVAP*(-1. ) 

ELSE 

anfilt-bncap 

PLUX(1)-ANFILT 

DETArN-DETAlN- (EVAP*AF ) 

END  IP 
C 

ENDIF 

C 

C 

C  CHANGES  IN  SOIL  MOISTURE  CONTENT 

C  . 

c 

c 

c 

SWP(NLL)  — 102.0 
DO  230  I-l.NL 

C  If  SWF  in  cell  is  greater  then  0,  it  is  saturated  and  flux  must 
C  therefore  be  0. 

IF(SWF(I+1) .GE.O.O)FLUX(I+1)-O.Q 

C  ANFLUX  represents  the  net  change  in  moisture  content  in  the  ceLL. 
ANFLl»X<  I  )-FLUX(  I )  -FLUX  ( 1  +  1 ) 

ANFLUX <1 ) -ANFLUX ( I ) *AF 

C  Recalculate  theta  according  to  the  change  influx(per  unit  area). 

THETA (I )-( VOL( I )+ANFLUX( I ) ) /TCOM< I ) 

C  Due  to  recalculation,  theta  may  be  greater  than  possible  water  content 

C  at  saturation  and  therefore  it  is  necessary  to  reset  SWP  to 

C  0  and  thata  to  the  water  content  at  saturation,  the  value  of  which  is 
C  entered  into  the  model. 

IF  <THETA(I).GE.SR1.AND.I.LE.NLA)SWP(I)-0.0 
IF  (THETA(I) .GE.SR2  AND. I .GT . NLa. AND . I .LE . NLH)SW?( I )-0 . 0 
IF(THETA(I).GE.SR3.AND.I.GT.NLH)SWP(I)-0.0 
IF(THETA( I ) . GE . SRI . AND . I . LE . NLA)THETA{ I )-SRl 
IF  f  THETA< I ) . GE . SR2 . AND , I . GT . NLA . AND . I . LE . NLH  >THETA( I >-SR2 
230  IF(THETA(I ) .GE . SR3 . AND . I . GT.NLB)THETA(I)-SR3 

C 
C 
C 

c 


£ - jr 
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C  CALCULATE  CUMULATIVE  TOTALS 

C  - - - 

c 

c 

c 

CUMDRN-CUMDRN+FLUX ( NLL ) *  AF 

evapi-evap*af+evapi 

CINFIL-CINFIL+ANFILT*AF 

S0G-THETACD/SR1 

C 

c 

c 

c 

c 

C 

c 

c 

c 

c 

C  . * . . . 

c  TERMINAL  SECTION  WRITE  OUT 

C  . 

c 

c 

c 

c 

c 

C  To  print  out  data  for  ovary  timo  incroaent  for  which  PPT  data  is 
C  entered,  check  ICOUNT  to  see  if  that  period  has  passed  by. 

IF (ICOUNT.LT. (DT*3600)>  GOTO  9995 
ICOUNT* 0 


C 

c 

c  CALCULATE  TIME  FROM  THE  START 

C  . 

T-T*AF/3600 
WRITE( 6 , 1170  )T 

1170  FORMAT ( ’ OSOIL  COLUMN  CONDITIONS  1 . F7 . 3 , IX. ’ HRS  SINCE 
&  SIMULATION  BEGAN’/) 

IF (TG.EQ. 06*00 . 0 )TG*0 . 0 


C 

C 

c  WRITE -OUT  CONDITIONS  OF  SOIL  COLUMN 

c  . . . . . . 

c 

c 

IF ( IOUT . EQ . 0 )GOTO  305 


WRITE(6 , 7700 ) 

7780  FORMAT ( ’  CeLl  Depth  SWP  Theta  Hyd  cond  Net’, IX. 
A'fluac  Rel  sat') 

DO  300  I-l.NL 

IF ( I . LE . NLA )SOG-THETA ( I ) /SRI 

IF ( I . GT . NLA . AND . I . LT . NLBB )SOG-THETA ( I ) /SR2 

IF (I . GE . NLBB )SOG-THETA( I ) /SR3 

300  WR1TE/6. 1190)1. DEPTH(I ) . SWP(I ) . THETAC I ) . COND( I) . ANFLUX ( I ) . SOG 

1190  FORMAT ( 16, 3F8 . * , 2F1* . 9 , F9 . 3 ) 

C 


JtL _ L 


c 

C  WATER  BALANCE  CHECK 

C  . 

c 

C  Philips  (1964)  simple  water  balance; 

C - 

C 

c 

C  Amount  added 

C  (Initial  soil) -(Current  soil)  -by  -  Evaporation- 

C  (  moisture  )  (  moisture  )  infiltration  loss 

C 

305  WATN-0. 

DO  310  I-l.NL 

310  WATN«TCCM(I)*THETA(I)+WATN 

B  AL-WATN -WAT  I  -CINFIL+EVAPI+CUMDRN 
WRIT£(6, 1200 )BAL 

1200  F0RMAT( ’ OBalance  check  on  soil  column  water  status  "'.F12 
8AL-(BAL*100. )/WATN 
WRITE (6, 1210 )BAL 

1210  FORMAT ( '  Balance  check  as  column  water  vol.  ■',F12.7,' 

C 

C 

IF ( IOUT . EQ . 0 )GOTO  306 

WRITE ( 6 , 1220 ) EVAPI . PPTT , CINFIL . CUMDRN 

1220  FORMAT ( *  Cumulative  evaporation  ■  1 ,F12.8/ 

&'  Cumulative  precipitation  •  ’,F8.4/ 

&'  Cumulative  infiltration  -  '.F10.6/ 

&’  Cumulative  drainage  »  ’ ,F10.6/) 

306  IF ( DETAIN. EQ.DETCAP) THEN 

WRITE(6, 1222) 

1222  FORMATC  Detention  capacity  exceeded’) 

WRITE ( 6 , 1230 )RT0T , RTOT/ . 0254 , T 

1230  FORMAT('  Runoff  total  in  the  last  period' ,F10. 7, 2X, 'm’ / 
&  '  Runoff  total  in  the  last  period’ ,F10 . 7 , 2X, ' ins ’ , 

$  F7 . 3/ ) 

ELSE 

WRITE(6, 122DDETAIN 

1221  FORMATC  Surface  water  -  \F10.6) 

WRITE(6, 1226) 

1226  FORMATC  No  runoff’) 

END  IF 
C 
C 
C 

C  CREATION  OF  ARRAY  DATA 

C  . . 

C 

c 

C  Runoff  is  recorded  in  array  WDATA 

C  The  runoff  for  each  soil  column  is  weighted  according  to  the 
C  percentage  area  which  it  occupies  in  the  catchment  area 
iww*^W 

WDATA ( . i www ) - ( RTOT / . 02 54 ) * ( I PC AREA / 100 . ) 

RTOT-O . 0 

■  ■  -  ■  «***  ■  i 

rrrmTrT  i 


Drainage 

loss 


.7) 


Z’/) 


9995  CONTINUE 
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C  End  of  simulation  of  single  soil  column,  it  mora  than  ona,  than  return  to 
C  to  the  beginning  of  thia  subroutine  to  repeat  for  next  soil  column 

34543  CONTINUE 

DO  76567  S-l.MM 

C  S*aa  the  weighted  runoff  for  each  soil  column  to  derive  total  runoff 
C  passed  back  to  CMPHYD  as  DATA 
CUMDATA-0 . 

DO  54345  J-l.NSCOL 

CUMDATA-WDATA (I , J) +CUMDATA 
54345  CONTINUE 

DATA(I)-CUMDATA 
76567  CONTINUE 


IR-f*W-l 

RETURN 

END 


SUBROUTINE  HYDCON(X , SATCON , SR, Z, Y) 

C  This  subroutine  calculates  hydraulic  conductivity  for  each  layer 
C  from  the  given  soil  moisture  characteristic  curve. 

C  Uses  the  Millington  and  Quirk  method 

DIMENSION  X(20).Y(20),Z(20) 

DO  845  1-1,20 

IIJ-20-I+1 

XII-X(IIJ) 

TOPS-Q . 

BOTS-O . 

DO  046  J-1,20 
JF-20- J+l 
YJJ-YCJF) 

046  BOTS-((2*J-l)*YJJ**(-2))+BOTS 

II-I 

DO  047  J-11,20 
JF-20" J+l 
YJJ-Y(JF) 

847  TOPS-C (2*J+1-2*I )*YJJ#*( -2) >+TOPS 

JT-20-I+1 

845  Z( JT)-SATC0N* (X( II ) /SR )*TOPS/BOTS 
RETURN 
END 


SUBROUT I NE  TWO ( NA , NB . THETA . X . SWP , Y , G , HPOT . DEPTH . GZ . COND . Z ) 

C  This  subroutine  calculates  soil  water  pressure,  hydraulic  potential 
C  and  hydraulic  conductivity  for  each  cell  as  conditions  change 
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C  during  simulation. 


DIMENSION  THETA (20) ,X(20) , SWP(20) ,¥(20) ,G(20) . HPOT(20) , 
&DEPTH(20),GZ(20),COND(20),Z(20) 

DO  15  I-NA.NB 
DO  16  J-1,19 

IF(THETA(I).GE.X(J).AND.THETA(I) .LT.X(J+1))SWP(I)-Y(J)+G(J)* 
i  (THETA(I)-X(J)) 

16  CONTINUE 

HPOT ( I ) -SWP  ( I ) -DEPTH ( I ) 

DO  17  J-1,19 

IF(THETA( I) .GT .X( J) .AND . THETA(I) . LE .X( J+l ) )COND(I )-Z( J)+GZ( J)* 
&  (THETA(I)-X(J)) 

17  CONTINUE 

15  CONTINUE 
RETURN 
END 


SUBROUTINE  GRAD(G,GZ, Y.X, Z) 


C  This  subroutine  calculates  the  gradients  of  the  suet  ion -moisture 
C  and  hydraulic  conductivity-moisture  curves. 

C 

DIMENSION  G(20) ,GZ(20) ,Y(20) ,X(20) ,Z(20) 

DO  261  1-1,19 

G(I)-(Y(I+1)-Y(I))/(X(I+1)-X(I)) 

261  GZ(I)-(Z(I+1)-Z(I))/(X(I+1)-X(I)) 

RETURN 

END 


SUBROUTINE  SMCURV ( SR , NQ , AX , Y , XNEW . YNEW . SCURV ) 


C  Generates  a  stochastic  suction  moisture  curve  to  be  fed  into 
C  soil  moisture  model 
c 
C 
C 

DOUBLE  PRECISION  G05DDF 
DOUBLE  PRECISION  AX, SCURV 

DIMENSION  AX( 20 ) ,X( 20 ) ,XNEW( 20 ) , YNEW (20 ),G(20),Y(20) 

C 

C 

C  Determine  the  stochastic  values  of  moisture 
C 

X ( 1 )-C0 5DDF (AX ( 1 ), SCURV ) 

IF (X( 1 ) . LT . 0 . )X( 1 )-0 . 001 
C 

DO  100  1-2. NQ 
X(I)^305DDF(AX(I) .SCURV) 

100  IF (X(I).LE.X(I-i))X(I)-X (I'l) +0.001 
IF(X(NQ) . GE. SR )SR-X(NQ) +0.001 


C 

C  Calculate  gradients  of  this  new  suction-moisture  curve 
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NNQ-NQ-l 
DO  200  I-l.NNO 

200  G(I)-(Y(I+1)-YCI>>/(XU+1)-X<I)) 

C 

C  Calculate  max  and  min  moistura  values,  and  determine  the  size  of 
C  equal  intervals. 

C 

XMAX-RMAX(X.NQ) 

XMIH-K  IH(X.NQ) 

XINT-  J^AX-XMIN)/19. 

C 

C  Determine  the  new  values  of  moisture-equal  intervals 
C 

XNEW(1)-XMIN 
DO  300  1-2,19 

300  XNEW(I  )-XNEW(  1)+(XINT*  (I~l) ) 

XNEW(20 )-XMAX 
C 

C  Determine  the  associated  new  values  of  suction 
C 

DO  350  1-1,19 
DO  400  J-l.NNQ 

IF(XNZW(I).GE.X(J).ANDXNEW(I).LT.X(J+1)) 

&  YNEW(I)-Y(J)-H3(J)*(XNEW(I)'X(J)) 

400  CONTINUE 
350  CONTINUE 

YNEW(20)-Y(NQ) 

C 

C 

C 

C 

RETURN 

END 


FUNCTION  RMAX  (X.NQ) 


C  Determines  the  maximum  real  in  an  array 


DIMENSION  X(NQ) 

C 

RMAX-XU) 

DO  10  1-2 ,NQ 

10  IF<X( I ) .GT . RMAX )RMAX-X( I ) 

C 

RETURN 

END 


FUNCTION  RMIN(X.NQ) 


C  Determines  minimum  real  in  an  array 


DIMENSION  X(N0) 
C 


RMIN-XC1) 


10 

c 


DO  10  1-2, HQ 

IF(X(I) .LT.RMIN)RMIN-X(I) 

RETURN 

END 


SUBROUTINE  PRTHYD 

C  CONVERTS  Q  HYDROGRAPH  TO  STAGE  HYDROGRAPH  FOR 
C  SPECIFIED  CROSS  SECTION 

C  ID-Q  BYD  INPUT 

C  IDR-CROSS  SECTION  ID 

C  NFK-2  OR  CREATOR  FOR  CONVERSION  Q/ STAGE 
C  NPK-1  Q  HYD 

C  NPK-0  0  PEAK  AND  VOLUME  ONLY 

C  IN  -  FORMAT  OF  OUTPUT 

C  IN  -0  REGULAR  FORMAT 

C  IN-1  PRINT  DISCHARGE  ONLY  IN  SINGLE  ENTRY  PER  LINE 
C  THIS  SUBROUTINE  PRINTS  THE  COORDINATES  OF  A  HYDROGRAPH. 

COMON/ BLOCKI/  OCFSOOO,  6)  ,DATA(311)  ,CFS(300)  ,CTBLE(50 , 11 )  , 
&RAIN(300).ROIN(6). 

&A(20 , 70 ) ,Q(20 , 70 ) , DEEP (20 ,  70 ) , ITBLEC  50 .2) , DP(20 ) , SCFS (20 ) , C (20 
&EALPHA(20) . IEND(6) ,DA(6) ,DIST(6) , SEGN ( 6 ) , DT ( 6 ) , PEAK ( 6 ) , ISG(6) , 
&NPU ,  NBD ,  NER .  MAXNO .  NO*W ,  ICC ,  NCODE ,  TIME ,  KCODE ,  ICODE 
COMON / BL0CK2 /  PERQ(20 , 70 )  ,TQ(20 . 6)  ,CC(20)  ,LL(6) .  INRC.LRC 
DIMENSION  DUM4Y(300) ,S(300,6), PEAKS 
DIMENSION  ISG(6) 

C  New  variables  used 

C  S  stage  equivalent  of  OCFS 

C  PEAKS  peak  stage  (equivalent  of  PEAK) 

C  Input  data  is  read  into  the  subroutine. 


ID-DATA(l) 

NPK-DATA(2) 

IDR-0ATA<3) 

IN-DATA U ) 

^IEND(ID) 

WRITE(6.40)ID.NPK 
TIME  1-0 

IFCNPK.LT.  DGOTO  32 
IF(NPK.LT.2)GOTO  2 
:  CONVERSION  TO  STAGE  HYDROGRAPH 

:  CHECK  RATING  CURVE  ENTERED 

IF(IDR.Eg.0)THEN 

WRITE ( 6 , * ) ' NEED  TO  ENTER  RATING  CURVE  ID' 

RETURN 

ENDIF 

C  CHECK  IF  MULTIPLE  ROUTING  INVOKED 

IF( IDR.GT . 6)GOTO  51 
IF (TQ (10, 1DR)  . GT.0)THEN 
DO  50  1-1,20 
50  0(I,IDR)-TQ(I. IDR) 

ENDIF 
JJ-IDR 
GOTO  7 


i 


f 

> 


t 


« 

i 


I 


j 

* 


v 


i 


]> 
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C  SEGMENT  HYDROGRAPH 

51  JJ-IDR/10 

7  DO  3  I-l.M 
J-l 

6  IF(OCFS(I, ID) .LE.Q(J, IDR) )GOTO  4 

J-J+l 

IF(J.GT.20)THEN 

WRITE ( 6, *> 'RATING  CURVE  EXCEEDED,  STOPPED’ 

RETURN 
END  IF 
GOTO  6 

4  IF(OCFS(I, ID) . EQ.Q(J, IDR) ) THEN 

S{X(X0)«J.JJ) 

GOTO  3 
ENDIF 

C  INTERPOLATE 

S(I ,  ID)KJ(  J-l,  JJ)+(  ( (OCFS(I , ID)-Q(J-1, IDR) )*(C(J,JJ)_ 
&C(J-1.JJ)))/(Q(J.IDR)-Q<J-1,IDR))) 

3  CONTINUE 

C  TIME  ARRAY 

2  DO  8  I-l.M 

DATAC I) -TIME 1 

8  TIME1-TIME1+DT ( ID ) 

J-0 

M4-M+4 

M5-W4/5 

IF(NPK . LT . 2)G0T0  27 
IF ( ICODE . EQ . 0 ) THEN 
WRITE(6 , 9) 

GOTO  10 
ENDIF 

WRITE( 6 , 11 ) 

10  IF(IN.GT.O)TREN 

IF ( ICODE. EQ.O) THEN 
DO  38  I-l.M 

38  WRITE ( 6 , 28 )S( I , ID) 

RETURN 

ENDIF 

DO  43  I-l.M 
SCI. ID)-S(I.ID)*0. 3048 
43  WRITEC6 , 28 )S( I , ID) 

RETURN 

ENDIF 

IF < ICODE. GT.O) THEN 
DO  45  I-l.M 

45  SCI, ID) -SCI, ID) *0.3048 

ENDIF 

39  J-J+l 

WRITE ( 6 , 30) (DATAC I ) , S( I , ID) , I-J ,M,M5) 

IFC J-M5 ) 39 , 13 , 13 
13  ROINl-ROINC ID) 

DO  16  1-1.20 

IFCQC I . IDR ) *PEAK (ID) ) 16 . 17 , 17 

16  CONTINUE 

17  IF(Q(I , IDR) . EQ. PEAK( ID) )THEN 
PEAKS<(I,  JJ) 

GOTO  18 
ENDIF 


I 

I 

i 

I 

f 


L 
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PEAKS-PEAK(ID)*((C(I-l,JJ)-C(IIJJ))/(Q(I-l,IDR) 

&-Q( I , IDR) ) )+C(I-l, IDR)*(C(I~1, JJ)- 
&C(I,JJ)))/(Q(I-1.IDR)-Q(I,IDR))) 

18  IF ( ICODE . EQ . 0 ) THEN 
WRITE ( 6 , 14 )ROIN 1 , PEAKS 
RETURN 

ENDIF 

PEAKS-PEAKS*0.3048 
RQIN1-ROIN(ID)*0. 0283168 
WRITE (6, 15 )ROINl, PEAKS 
RETURN 

C  DISCHARGE  HYDROGRAPHS 

27  IF ( ICODE. EQ.l) THEN 

C  METRIC 

WRITE(6, 21) 

DO  23  I-l.M 

23  Dt*fff(  I  )-OCFS(  I,  ID  >*0.0283168 

PEAK1-PEAK(ID)*0 . 0283160 
ROINl-ROIN( TD  >  *0 . oZBolSS 
GOTO  20 
ENDIF 

C  IMPERIAL 

19  WRITE( 6 , 25) 

DO  26  I-l.M 

26  DUl*«CI)-OCFS(I.ID) 

PEAK  1- PEAK  ( ID ) 

ROINI-ROIN(ID) 

20  IF(IN.GT.0)THEN 
DO  29  I-l.M 

29  WRITE(6t28)DUNMY(I) 

RETURN 

ENDIF 

31  J-J+l 

WRITE(6.30)(DATA(I),DUNWY(I).I-J.M.M5) 

IF( J-M5)31, 32, 32 

32  IF ( ICODE . ME . 0 )GOTO  34 
ROINl-ROIN(ID) 

PEAK 1-PEAK (ID) 

WRITE(6, 35)R0IN1 , PEAK1 
RETURN 

34  ROIN1-ROIN( ID )*0 . 0283168 
PEAK1-PEAK( ID ) *0 . 0283 168 
WRITE ( 6 , 36 )ROINl . PEAK1 
RETURN 

21  FORMAT 1 10X,  "TIME" ,  6X,  "  FLOW” , 11X, "TIME” . 6X, "  FLOW ' , 11X, "TIME" . 
&6X , "FLOW" . 1 IX . "TIME" . 6X . "FLOW" , 11X , "TIME" . 6X .  "FLCW" /I IX . "HRS" . 
&7X , "  MS" , 12X , "HRS" ,  7X .  "  MS" . 12X. "HRS" . 7X. "  MS" . 12X . "HRS " . 

&7X , "  MS”. 12X."HRS".7X,"  MS") 

25  FORMAT (10X." TIME ",6X."  FLOW" , 1 IX. "TIME" , 6X. "  FLOW" . 11X . "TIME  . 

&6X . "FLOW" . 11X . "TIME" . 6X . "FLCW" . 11X , "TIME" . «X . "FLCW" / 1 IX . "HRS" . 
&7X , ”  CFS  " , 10X , "HRS" , 7X , "  CFS  " . 10X, "HRS" . 7X . "  CFS  " , 1CX . "HRS " . 
&7X ,  '*  CFS  "  ,  10X ,  "HRS"  .  7X,  "  CFS  ") 

30  FORMAT  ( 5< 5X, F10 . 3 , F10 . 3 ) ) 

40  FORMAT (  ’PRINT  HYD1  .T21,  1  ID-MI. T29, 'NPK-' .11) 

36  FORMAT( 1H0.9X. "HYDROGRAPH  VOLUME-", F20.0."  CUMEC  "/10X.PEAK 
&  DISCHARGE  RATE  -".FIO.O, "CMS”///) 

35  FCRMAT(1H0,9X, " HYDROGRAPH  VOLUME-" . F20 . 0 . "  CF  "/I OX, "PEAK 
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14 

15 

11 

9 

) 

9 

28 


4  DISCHARGE  RATE-" ,F10 . 0 , "CFS"///) 

FORMAT { 1H0 , 9X,  "HYDROGRAFB  VOLUME-"  ,  F20 . 0,  "  CF  "/10X,"PEAK 
&  ELEVATION  -".F10.0."  FEET"///) 

FORMATUH0.9X,  "HYDROGRAPH  VOLUME-" ,  F20 .0  ,  "CUMECS"/10X.  "PEAK 
&  ELEVATION  -" ,F10. 0 , "METRES"/// > 

FORMAT! 10X ,  "TIME" ,  SX ,  "ELEV”  ,  UX ,  "TIME" .  6X.  "ELEV" .  UX ,  "TIME” , 

St  6X.  "ELEV",  UX,  "TIME",  6X,  "ELEV",  I1X,  "TIME”.  6X,  "ELEV", /11X,  "HRS". 

St  7X ,  "M  " .  12X,  "HRS"  ,  7X,"M  " ,  12X,  "HRS"  ,  7X.  "M  ",I2X,"HRS", 

&  7X,"M  " ,  12X,  "HRS"  ,  7X,  "M  ") 

FORMAT { 10X,  "TIME" , 6X ,  "ELEV"  ,  UX ,  "TIME" ,  6X ,  "ELEV" .  UX.  "TIME" . 

St  SX ,  "ELEV"  ,  UX ,  "TIME"  ,  SX ,  "ELEV"  ,  11X ,  "TIME"  , 6X,  "ELEV"  ,  /UX .  "HRS" . 
St  7X,  "FT" .  12X ,  "HRS"  ,  7X.  "FT"  ,  12X,  "HRS"  ,  7X,  "FT”  ,  12X,  "HRS"  , 

St  7X,  "FT" ,  12X,  "HRS" ,  7X,  "FT") 

FORMAT  <F  10. 3) 

END 


SUBROUTINE  PUHYD 

C  THIS  SUBROUTINE  PUNCHES  HYDROGRAPHS  IN  FORM  TO  BE  USED  BY 
C  SUBROUTINE  RECHD 

CCXKON/HLOCKl/  OCFS(300,6) ,DATA(311) ,CFS(300) ,CTBLE(50, 11). 
4RAINC300) ,ROIN<6) , 

4AC20 , 70)  ,0(20 , 70) ,DE£P<20 , 70) . ITBLE( 50 ,2) ,DP(20) , SCFSC20) , C 1 20 ,67, 
RZALPHAC20) , IEND(6) ,DA(6) ,DIST(6) , SEGN(6) ,DT(6) ,PEAK(6> , ISG(6 ) , 

RNPU, NHD . NER , MAXNO , NCOfK,  ICC , NCODE , TIME . KCODE ,  ICODE 
DIMENSION  OUfKYOOO) 

1D-DATM1) 

M-IEND(ID) 

IF < ICODE. EQ.OiGO  TO  3 
DA1-DA(ID)*2.590 
PEAK1*PEAK( ID )*0 . 02832 
ROIN1-ROIN ( ID ) *0 . 02832 
DO  *  I-l.M 

DUSKY <  I ) -OCFS  <  I .  I D )  * 0 . 02 832 

4  CONTINUE 

WRITE! 7 , 5IID.NHD, DTCID) .DAI, PEAK1 .ROIN1 , 1  END ( ID) , ICODE 
WRITE(  7 . 6 )  ( DUSKY  (I), I-l.M) 

RETURN 

3  WRITE ( 7 , 1 )ID.NHD . DTCID) ,DA( ID) , PEAK ( ID) ,ROIN(ID) . IEND( ID) , ICODE 
WRITE  (7.2)  (OCFS(I. ID), I-l.M) 

RETURN 

C 

1  FORMAK  'RECALL  HYD'  ,T21,  'ID-'  ,11.  T29,  'HYD  NO’,I),Tl2,  'DT-'  ,F9. 

AS,'  HRS’ ,T61. 'DA-' ,F8  3, '  SQ  MI '/Til , ' PEAK-' , F7 . 0 . 'CFS' ,T40 . 'RO- ' , 
AF20.0,"  CF  ” . T60 , "NO  PTS-" . I3/21X , "CODE-" . U/T21 . 

4'FLOW  RATES") 

5  FORMAT!  'RECALL  HYD '  ,T2I .  '  ID-'  .  II. T29 .  ' HYD  NO'  ,  13 , T42 .  ' DT-'  .  F9 

46,’  HRS'  ,T61.  'DA-'  ,F8.3. '  SQ  KM7T21. 'PEAK-' ,F7.2. 'CMS' ,T40. 'RO-' . 
4F20 . 0 , "  CUMEC  ".TSO.  'NO  PTS-" . I3/21X. "CODE-" . I1/T21 . 

4  FLOW  RATES') 

2  FORMAT  (T21.7F8.0) 

6  FORMAT  (T21.7F8.2) 

END 


A 
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SUBROUTINE  HPLOT 

C  THIS  SUBROUTINE  PLOTS  EITHER  1  OR  2  HYDROGRAPBS  ON  A  SET  OF  AXIS 

COMON  /  BLOCK  1  /  OCFS  ( 3  00 , 6 ) , D ATA ( 3 1 1 ) . CFS ( 300 ) ,  CTBLE  (50,11), 

&RAIN ( 300) ,R0IN(6) , 

&A(20 , 70 ) ,0(20 , 70 ) , DEEP (20 . 70) , ITBLE( 50,2) , DP(20 ) , SCFS (20 ) , C(20 . 6 ) , 
&ZALPHAC20) ,IEND(6) ,DA(6),DIST(6),SEGN(6),DT(6) ,PEAK(6),ISG<6), 

&NPU , NHD . NER , MAXNO . NCCMM , ICC , NCODE .TIME, KCODE . ICODE 
IDl-DATA(l) 

ID2»DATA(2) 

DATA  ZERO.  PLUS.  BLANK,  DASH,  DOT/ ' 0 

MAX-121 

J-l 

C  ARE  THERE  1  OR  2  HYDROGRAPHS 

IF  (ID2)  1,1.2 

C  DETERMINE  HIGHEST  PEAK  IF  2  HYDROGRAPHS 

1  QMAX-PEAK(IDl) 

GO  TO  14 

2  IF  ( PEAK( ID! ) -PEAK ( ID2 > )  3,3.4 

3  QMAX-PEAKC ID2) 

GO  TO  5 

4  QMAX-PEAKC ID  1) 

C  IF  2  HYDROGRAPHS  DETERMINE  LARGEST  DT  AND  INTERPOLATE  OTHER 

C  HYDROGRAPH  IF  NECESSARY 

5  IF  (DT(ID1)-DT(ID2) )  6,13,7 

6  L-ID1 
K-ID2 
GO  TO  8 

7  L-ID2 
K-ID1 

8  M-IEND(L) 

TID-DT(K) 

TIDH-0. 

DO  11  1-2, M 
TIDH-TIDH+DT(L> 

IF  (TID-TIDH)  10.9.11 

9  J-J+l 
CFS(J)-OCFS(I.L) 

TID-TID+DKK) 

GO  TO  11 

10  J-J+l 

CFS(J)-OCFS(I-l.L)^((TID-TIDH*DT(L)  )/DT(L)  ) * (OCFSC I . L ) -OCFS ( I  - i . L 1 
&) 

TID-TID+DKK) 

11  CONTINUE 
IEND(L)-J 
DT(L)-DT(K) 

DO  12  1-2. J 

12  OCFSU.L)-CFS(I) 

13  IF  ( I END ( IDl ) - IENDC ID2) ) 

U  IENDC IDl) 

GO  TO  16 

15  M-rEND(ID2) 

IE  XM  -  M 

C  DETERMINE  TIME  3CALF 

XCCL  -  XM  121 
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YSCL^MAX/50. 
c  PLOT  HYDROGRAPHS 

DO  20  I-l.MAX 

20  CFS(I)«DASH 

IT ( ICODE. EQ.O)GO  TO  *9 
WRITE<6,50) 

50  FORMAT (T2, "FLOW  RATE  (CMS)") 

QMAX1-CMAX*0 . 02832 

WRITE ( 6 . 4 1 )QMAX 1 , DOT , (CFS ( I ) , 1-1 . MAX ) , DOT 
GO  TO  51 

49  WRITE(6,40) 

48  FORMAT ( T2, ’FLCW  RATE  (CFS>’> 

WRITE (6,41 )QMAX . DOT , (CFS(I) . I-l.MAX) .DOT 

51  Q1KJMAX 
Jl-10 

DO  37  J-l , 50 
IF  (J-Jl)  23,21.23 

21  DO  22  I-l.MAX 

22  CFS( I )-DASH 
GO  TO  25 

23  DO  24  I-l.MAX 

24  CFS(I)-BLANK 

25  Q2-Q1-YSCL 
DO  28  1-2, M 

IF  (OCFS(I.TDl)-Ql)  26,27.28 

26  IF  (OCFS(I,IDl)'Q2)  28.28.27 

27  XI  -  I 

K  -  XI  /  XSCL  +  1. 

CFS(K)-ZERO 

28  CONTINUE 

WRITE  (6.44)  DOT, (CFS( I). I-l.MAX). DOT 
IF  (ID2)  34.34.29 

29  DO  18  I  -  1.  MAX 
18  CFS(I)  -  BLANK 

DO  33  I-l.M 

IF  (0CFS(I.ID2)‘Q1)  30.3 1,33 

30  IF  (OCFS ( I , ID2 ) *02 )  33.33.31 

31  XI  -  I 

K  -  XI  /  XSCL  1. 

CFS(K)-PLUS 

33  CONTINUE 

WRITE  (6.42)  (CFS(I>. I-l.MAX) 

34  IF  (J-Jl)  36.J5.36 

35  Jl-Jl+10 

IF ( ICODE. EQ.C)GO  TO  52 
QDH32*0. 02832 
WRITE(6,43)0D 
GO  TO  36 

52  WRITE(6, 43)02 

36  Q1-Q2 

37  CONTINUE 
CFS ( D-TIME 

DTT-DT ( ID  1 )* <  XM  -  1. >  /  12. 

C  PUT  TIME  ARRAY  IN  CFS  AND  WRITE  TIME  SCALE 
DO  38  1-2.13 

35  CFS( I  )-CFS( I - 1  )+D7T 

WRITE  (  6  .  <•  5 )  fCFSd). 1-1.13) 

WRITE 
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IF  (NPU)  40,40.39 

39  WRITE  (7,47)  ID1.ID2 

40  RETURN 
C 

41  FORMAT ( 1X,F7 . 0 , 123A1) 

42  FORMAT ( 1H+ , 8X . 12 1A1 ) 

43  FORMAT  (1H+.F7.0) 

44  FORMAT (8X.123A1) 

45  FORMAT(T3 , 13F10 .2) 

46  FORMAT ( 4 9X, 'TIME  HOURS’///) 

47  FORMAT (  ’PLOT  BYD*  ,T21,  ’  ID  I-\ II , T29, ' ID  II-' ,11) 

END 


SUBROUTINE  adhyd 

C  THIS  SUBROUTINE  ADDS  TWO  HYDROGRAPHS. 

C 

C  INCLUDES  CORRECTION 

C  IF  DT(IDl)  AND  DT(ID2)  ARE  DIFFERENT  "AND" 

C  ID  IS  EITHER  EQUAL  TO  ID1  OR  ID2  THEN 

C  EXSISTING  CODE  WOULD  HAVE  LOST  ITERATION  PERIOD  OF  ONE  OF 
C  THE  INFLOW  HYDS  -  THIS  IS  NOW  CORRECTED 

CCmDN/ BLOCK  1/  OCFS(300 . 6) ,  DATA(311)  ,CFS(300 )  ,CTBL£(  50 , 11 ) , 
&RAIN(300 ) ,ROIN(6 ) , 

&A(20 , 70 ) , Q(20 , 70) .DEEP (20 , 70 ) , ITBLEC  50 , 2 ) , DP(20 ) . SCFS (20 ) , C{ 20 , 6 ) , 
&ZALPHA( 20 ) , IEND( 6 ),DA(6) ,DIST(6) .SEGNC6 ) , DT( 6 ) , PEAK ( 6 ) , ISG(6) . 
&NPU ,  NHD ,  NER ,  MAXNO ,  NCOM,  ICC ,  NCODE ,  TIME  ,  KCODE .  ICODE 

ID-DATA(l) 

NHD-DATA(2) 

IDl-DATAO) 

ID2-DATA(4) 

KK-0 

C  CHECK  ARRAYS  ARE  NOT  EMPTY 

IF ( IEND{ ID1 ) . EQ. 0 .OR . I END ( ID2) . EQ. 0 ) THEN 
WRITE(6.*) 'ONE  HYDROGRAPH  BEING  ADDED  IS  ZERO' 

IF  (TEND  (TDD  EQ.0)THEN 

K-ID2 

GOTO  53 

END  IF 

K*ID1 

53  DO  52  I-l.IEND(K) 

52  OCFS(I.ID)-OCFS(I,K) 

PEAK( ID)-PEAK(K) 

ROIN( ID)-ROIN(K) 

DA(ID)-DA(K) 

IEND(ID)-IEND(K) 

DT(ID)»DT(K) 

GOTO  27 
ENDIF 

IF(DTdDl)  .EQ.DT(ID2))GOTO  54 
IF ( I D  NE  IDI  AND . ID  NE  ID2 )GOTO  54 
C  DANGER  OF  CONFUSION  IN  DT .  ALTER  ID  TO  KK 
55  DO  56  KX-1.5 


a 


T 
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IF(KX.EQ. ID 1) GOTO  56 
IF(KK.EQ.ID2)GOTO  56 
GOTO  57 

56  CONTINUE 

57  ID-KK 

54  PEAK(ID)  -  1. 

C  MAKE  TIME  INCREMENTS  EQUAL  IF  NOT  EQUAL.  USE  SMALLER  INCREMENT 

IF  (DT(ID1)-DT(ID2))  1,3,2 

1  DTUD)-DT(IDl) 

L-ID1 

K-ID2 
GO  TO  6 

2  DT(ID)-DT(ID2) 

L-ID2 

K-ID1 
GO  TO  6 

3  DT(ID)-DT(ID1) 

IF  ( IEND ( ID1) - IEND ( ID2 ) )  4.4,5 

4  M3-IENDUD1) 

K1-ID2 

IEND(ID)-IEND(ID2) 

GO  TO  18 

5  M3-IEND(ID2) 

K1-ID1 

IEND(ID)-IEND(ID1) 

GO  TO  18 

C  DETERMINE  DURATIONS  OF  FLOW 

6  XIEND1-IEND(ID1)-1 
XIEND2-IEND(ID2)-1 
DUR1-XIEND1*DT(ID1) 

0UR2-XIEND2*DT ( ID2 ) 

IF  (DUR1-DUR2)  7,8,8 

7  IEND( ID)-DUR2/DT( ID )+l . 

M3-DUR1/DT(ID)+1. 

K1-ID2 

GO  TO  9 

8  IEND( ID )-DURl/DT ( ID)+1 . 

M3-DUR2/DT(ID)+1. 

K1-ID1 

9  IF  (IEND ( ID )-300 >  11,11,10 

10  IEND(ID)-300 

11  M2-IEND(K) 

J-l 

C  INTERPOLATE  ONE  HYDROGRAPH  IF  NECESSARY 

TIDH-0. 

TID-DT(ID) 

DO  15  >2, M2 
TIDH-TIDH+DT(K> 

12  IF  (TIDH-TID)  15. 13,14 

13  J-J+l 

DATA  (J)-OCFS(I.K) 

TID-T1D+DTCD) 

IF  (J-300)  15.16.16 

14  J-J+l 

DATA  (J)-OCFS(I-l.K}*(iriD-TIDH+DT(K))/CT(KJ  >*(OCFS(I.K)-OCFSrl-l. 
&K)) 

TID-TID^DT ( ID } 

IF  (J-300)  12. 16. 16 


15  CONTINUE 

16  IEND<K)-J 
DO  17  1-2, J 

17  OCFS ( I , K ) -DATA ( I ) 

18  MhlEND(ID) 

RO  -  0. 

C  ADD  HYDROGRAPHS 

C  CONVERT  KK  TO  ID 

IF(KK.GT.O)THEN 
ID-DATA(l) 

DT(ID)-DT(L) 

ENDIF 

DO  20  1*1, M3 

OCFS  ( I .  ID  XCFS  ( I ,  ID  1 ) +OCFS  ( I ,  ID2  ) 

IF  (OCFS(I , ID)  -  PEAK(ID) )  20,20,19 

19  PEAKdD)  -  OCFS  (I ,  ID ) 

20  RO  -  RO  +  OCFS (I, ID) 

DA( ID)-DA(ID1 )+DA( ID2 ) 

IF  (PEAK(ID)  -  PEAK(Kl) )  21,22,22 

21  PEAKdD)  -  PEAK(Kl) 

22  IF  (M-M3)  25,25.23 

23  M3  -  M3  +  1 

DO  24  I  -  M3. M 
OCFS(I.ID)  -  OCFS  (I ,K1) 

24  RO  -  RO  +  OCFS(I.ID) 

25  ROIN(ID)  -RO  *  DT (ID) *3600 
IF  (NPU)  27,27.26 

26  WRITE  (7,20)  ID.NHD, ID1, ID2 

27  RETURN 
C 

28  FORMAT (  ’ ADD  HYD ' . T21 . ' ID-'.Il . T29 . '  HYD  NO- • . 13 , T45 . ' ID  I- '.1 1 

&T60,  ’ID  II-M1) 

END 


SUBROUTINE  SRC 

C  THIS  SUBROUTINE  STORES  AN  ELEVATION  -  END  AREA  -  FLOW  TABLE. 

CCfWON / BLOCK 1 /  OCFS (300, 6) ,DATA(311 ) ,CFS(300) .CTBL£(50. 11 ) . 
&RAIN(300) ,ROIN(6) , 

&A( 20 , 70 ) , Q ( 20 , 70 ) , DEEP (20,70), ITBLE ( 50 , 2 ) , DP( 20 ) . SCFS ( 20 ) , C  <  20 . 6 ) 
&ZALPHA(20) , IEND( 6 ),DA(6).DIST(6), SEGN(6) , DT (6 ) , PEAK( 6) . ISG(6) . 
&NPU .  NHD .  NER ,  MAXNO .  NCOW ,  ICC ,  NCODE .  TIME .  KCODE  ,  ICODE 
ID-DATA(l) 

VS-DATA( 2 ) 

C  VALLEY  SECTION  NUMBER 

C  REMAINING  DATA  ARE  ELEVATION.  AREA,  AND  FLOW  FOR  EACH  POINT  OF 

C  THE  RATING  CURVE 

IF (KCODE. EQ.OJGO  TO  2 
J-3 

DO  3  1-1.20 

DATA( J )-DATA( J ) / 0 . 3048 
DATA( J*1 )-DATA( J+l )/0  093 
DATA ( J+2 ) -DATA ( J+2 ) /0  02832 
J-J+3 

3  CONTINUE 


2  EMIN-DATA(3) 

J-3 

DO  1  1-1,20 

DEEFCI, ID)-DATA( J)-EMIH 
A( I , ID)-DATA( J+l ) 

Q(I , ID)-DATA(J+2) 

J-J+3 

1  CONTINUE 

RETURN 
END 

SUBROUTINE  CMPRC 

C  THIS  SUBROUTINE  COMPUTES  THE  DISCHARGE  END-AREA  ELEVATION 

C  RELATIONSHIP  FOR  A  VAINLY  SECTION. 

C 

C  IF  MUTIPLE  ROUTING  INVOKED  - 

C  COMPUTES  SEPARATE  RATING  CURVES  FOR  EACH  SEGMENT 

C  ALSO  COMPUTES  Z  FLOW  AT  EACH  ELEVATION  FOR  SEPARATE  SEGMENTS 

C 

C  IF  TURBULENT  EXCHANGE  INVOKED 

C  COMPUTES  THE  RATING  CURVE  USING  REDEFINED  AREA  AND  WETTED 
C  PERIMETER  CALCULATION  -  KNIGHT  TECHNIQUE 

C  FOUR  OPTIONS 

C 

C  NOTE  — -  TURBULENT  EXCHANGE  REDEFINITIONS  USED  "ONLY  ,f 
C  FOR  OUT -OF -BANK  ELEVATIONS 

C 

C  MULTIPLE  ROUTING  AND  TURBULENT  EXCHANGE  OPERATES  INDEPENDANTLY 

COWON / BLOCK  1  /  OCFSOOO  .  6) ,  DATA<3X1)  ,CFS< 300 ) , CTBLEC 50 , 11 ) . 
&RAINOOO),ROIN(6), 

&A(20 , 70 ) ,Q(20 , 70) , DEEP (20 , 70 ) , ITBLE ( 50 , 2) , DP( 20 ) , SCFS (20 ) ,C ( 20 , 6 ) . 
&ZALPHA(20) ,IEND(6),DA(6).DIST(6) ,SEGN(6).DT(6) .PEAKC6) ,ISG(6), 

&NPU .  NHD ,  NER , MAXNO ,  NCOW.  ICC ,  NCOQE .  TIME ,  KCODE .  ICODE 
DIKLITT ION  f*W(6),W(6),XM(70) 

COWON / BL0CK2 /  PERQ(20 , 70) , TQ<20 , 6) , CC(20 ) ,LL (6 )  .  INRC.LRC 

C  new  variables  used 

C  fM*  used  to  mark  segments  for  turbulent  exchange 
C  l-le£t  floodplain  2-channel  3-right  floodplain 

C  W  width  of  channel  segment 

C  H  channel  depth 

C  XM  min  elev  in  segment 

C 

ID-DATA(l) 

C  STORAGE  LOCATION  NUMBER.  (1-6) 

IT-DATA<2) 

C  TURBULENT  EXCHANGE  INCLUSION 

MR-DATA(3) 

C  MULTIPLE  ROUTING  INCLUSION 

VS-0ATAU) 

C  VALLEY  SECTION  IDENTIFICATION  NUMBER. 

NSEG-0ATA<5> 

C  NUMBER  OF  SEGMENTS  IN  THE  VALLEY  SECTION 
IF ( KCODE. EQ. 0)GO  TO  26 
DATA(6)-DATA(6)/0  3048 
DATA(7)-DATA(7)/0  3048 
26  ELO-DATAC6) 
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EMAX-DATA ( 7 ) 

C  MAXIMUM  ELEVATION  FOR  COMPUTATIONS . 

SLOPE 1-DATA ( 8 ) 

C  CHANNEL  SLOPE . 

SLOPE2-DATA(9> 

C  FLOODPLAIN  SLOPE. 

DIF-(EMAX-ELO)/19. 

C(l.ZD)-£LO 
DO  1  1-2,20 

1  C(I.ID)-C((I-1).ID)+DIF 

C  SET  AREA  AND  DISCHARGE  ARRAYS  -  0. 

IF(MR.GT.O)  GOTO  53 
DO  2  1-1,20 
A(I,XD)*0. 

2  Q< I , ID)-0 . 

DO  76  K-l.NSEG 
W(K)-0 

76  f**!(K)-0 

GOTO  54 

53  DO  55  J-10*ID+1, 10*ID+NSEG 
DO  56  1-1,20 

A( I , J)-0 
QCI,J)-0 
TQCI , ID)-0 
56  PERQ ( I , J ) -0 

55  CONTINUE 

54  J-10 

C  READ  N  VALUES  AND  SEGMENT  BORDER  POINTS. 

DO  3  I-l.NSEG 
SEGN (I) -DATA (J) 

IF (KCODE . NE . 0 ) DATA ( J+l ) -DATA (J+l ) /0 . 3048 
DIST ( I )-DATA( J+l ) 

3  J-J+2 

C  REMAINING  DATA  ITEMS  ARE  DISTANCES  AND  ELEVATIONS. 
IFCKCODE . EQ . 0 )GO  TO  27 
DO  28  I-J.310 
DATA( I )-DATA( I ) /0 . 3048 
28  CONTINUE 

27  JJJ-J 

DO  6  I-l.NSEG 

4  J-J+2 

IF  (DATA(J)  -  DISK  I) )  4.5.5 

5  ISG(I)  -J+l 

6  CONTINUE 

C  COMPUTE  CHANNEL  WIDTH 

IFdT.LT.  1)GOTO  78 

J-10 

35  DO  61  K-l.NSEG 

SELEV-0 

IF(SEGN(K ) ) 62 , 63 , 63 
63  IF(K.E0.  DGOTO  65 

IF(SEGN(K~ 1)) 64, 65, 65 

65  IF (K . EQ . NSEG)GOTO  61 

IF(SEGN(K+1) )66 . 61 , 61 

66  IF(K.LT. 21GOTO  67 

C  LEFT  HAND  FLOODPLAIN 

GOTO  61 
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C  CHANNEL 

62  IF(K.EQ. 1. OR. K. EQ. NSEG)THEN 

WRITE(6 , 70) 

RETURN 

ENDIF 

W(K)-(DATA(ISG(K)-l)-DATA(ISG(K-l)-l))/2 

t**!(K)-2 

GOTO  61 

C  RIGHT  HAND  FLOODPLAIN 

64  NMM(K)-3 

GOTO  61 

C  LEFT  HAND  FLOODPLAIN 

67  M*1(K)-1 

61  CONTINUE 

C  COMPUTE  DISCHARGES  AND  END  AREAS  FOR  EACH  SEGMENT. 

78  IFCIT.NE. 0)THEN 

WRITE (6, 170  JIT 
ENDIF 

DO  22  K-I.NSEG 
J-JJJ 
JJJ1-JJJ+1 
IF  (SEGN(K)J  7,7,8 

7  SLOPE-SLOPE 1 
GO  TO  9 

8  SLOPE-SLOPE 2 

9  SLPN-1 . 486+SLOPE+* . 5 

C  COMPUTE  AREA  AND  DISCHARGE  FOR  SEGMENT. 

DO  21  1-2,20 
AA-0.' 

P-0. 

J-JJJ- 1 
DEP2-0 . 

10  J-J+2 

IF  (J-ISG(KJ)  12,12,11 

11  IF(AA-. 001)21, 21, 20 

12  IF(DATA( J) -C( I , ID) )  13,10,10 

13  DEP1-C(I,ID)-DATA(J) 

IF  (J-JJJ1)  16.16,14 

14  XL-DATA(J-l)-DATA<J-3) 

DEP3-ABS ( DATA ( J- 2 ) -DATA ( J) ) 

XL-XL+DEP1/DEP3 

15  AA-AA+XL+(DEPl+DEP2)/2. 

P-P+SQRT ( ( DEP1 -DEP2 ) **2+XL**2 ) 

16  DEP2-DEP1 
J-J+2 

IF  C J-ISG(K) )  17,17,20 

17  IF  (DATA( J) -C( I , ID ) )  18.18,19 

is  depik:(i,id)-data(J) 

XL-DATA(J-l)-DATA(J-3) 

GOTO  15 
19  DEP1-0 

XL-DATA ( J- 1 ) - DATA ( J -  3 ) 

DEP3-ABS(DATA( J-2)-DATA( J) ) 

XL-XL+DEP2/DEP3 

AA-AA+XL*(DEPl+DEP2)/2. 

P-P+SQRT { (DEP1-DEP2 ) **2+XL**2 ) 

DEP2-0 . 

GO  TO  10 
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C  CHECK  IF  TURBULENT  EXCHANGE  INVOKED 

C  CHECK  IF  OUT-OF-BANK 

20  IF(ff«(K).LT.  l)GOTO  98 
IF(h*W<K).EQ.l)GOTO  90 
IF(»«(K).EQ.3)GOTO  91 

C  CHANNEL 

C  CHECK  OUT-OF-BANK 

IF(C(I,ID).LE.DATA(ISG(K)).AND.C(I, 

&ID).LE.DATA(ISG(K-1)))GOTO  98 
B«(DATA(ISG(K))+DATA(ISG(K-l)))/2-ELO 
IF(IT.LL.2)G0T0  92 
C  AREA  METHOD  3  AND  * 

AA-AA/2+(W(K)*H) 

92  IF(IT.EQ. 1. OR. IT. EQ. 3) THEN 
C  WETTED  PERIMETER  METHOD  1  AND  3 

P-P- (2*(C(I. ID)-CC (1-1) . ID) ) )+2*H 
ENDIF 

IF(IT.EQ. * )THEN 

C  WETTED  PERIMETER  METHOD* 

P*P+(2*( (C(I , ID)-C( CI-1) , ID) )**2+W(K)**2)**0 . 5) 

ENDIF 
GOTO  98 

C  LEFT  HAND  FLOODPLAIN 

90  L-K+l 
GOTO  95 

C  RIGHT  HAND  FLOODPLAIN 

91  L-K-l 

95  IF (IT . LT . 3 )GOTO  96 
AA-AA+((C(I,ID)-C(d-l).ID))*W(L)/2) 

96  1F(IT.EQ.2)THEN 
P*P+(C(I ,  ID)-C(d-l).ID)) 

ENDIF 

98  R-AA/P 

IF(SEGNfv"  i.T .  0 . 0 )  THEN 
SGN— S  K) 

GOTO  130 
ENDIF 

SGN-SEGN(K) 

130  IF(MR.LT.l)  GOTO  37 

C  COMPUTE  SEPARATE  R. CURVES  FOR  EACH  SEGMENT 

II"10*ID+K 

Q( I . II )-Q( I , II )+AA*R** . 6667*SLPN/SGN 
Ad ,  II  )"A(  X ,  II  )+AA 
GOTO  21 

C  ADD  DISCHARGES  AND  AREAS  FOR  ALL  SEGMENTS  TO  OBTAIN  TOTALS  FOR 
C  VALLEY  SECTION. 

37  0(1 . ID)-Q<I, ID)+AA*R*V66667*SLPN/SGN 

A(I , ID)“A(I , ID )+AA 

21  CONTINUE 
JJJ-J-3 

22  CONTINUE 

IF( ICODE . EQ . 0 )GO  TO  29 
IF(MR.LT. l)GOTO  *0 
C  FIND  MIN  ELEV  IN  EACH  SEGMENT 

J-13+2*NS£G 
DO  80  ^l.NSEG 
IF(M.EQ.  DTHEN 
XM(10*ID+M)-DATA(ISG(M)) 
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GOTO  81 
ENDIF 

XM<10#ID+M)"DATA(ISG{M“1)  > 

81  IF(J.GT.ISG(M))GOTO  80 

IF ( DATA ( J) . LT . XM ( 10* ID+M ) ) THEN 
XM( 10*ID+M)-DATA(J> 

ENDIF 
J-J+2 
GOTO  81 
80  CONTINUE 

DO  *1  J-10*ID+1,10*ID+NSEG 

WRITE (6, 47 )J 

DO  42  1-1,20 

CK(I,ID>*0.3048 

Al-A< I , J)*0 . 093 

Q1-Q(I, J)*0. 02832 

DEEPd.JJ-Cd.IO-XM^) 

IF(DEEP(I, J) .LT.0)THEN 

DEEPd.J>“0 

ENDIF 

WRITE (6, 32)  C1.A1.Q1 
42  CONTINUE 
41  CONTINUE 
GOTO  49 

40  WRITE (6 , 31 )VS 
DO  30  1-1,20 
C1^CI,ID)*0.3048 
A1-A(I , ID)*0 . 093 
Q1-Q(I . ID)*0 . 02832 
DEEPd.ID)-C(I.ID)-ELO 
WRITE(6, 32)C1, Al.Ql 
30  CONTINUE 
RETURN 

29  IF(MR.LT. l)GOTO  43 

C  FIND  MIN  ElEV  IN  SEGMENT 

J-13+2*NSEG 
DO  82  ^l.NSEG 

IF(M.EQ.1)THEN 

XM( 10*ID+M)-DATA(ISG(M) ) 

GOTO  83 
ENDIF 

XM( 10*ID+M)-DATA( ISG(M- 1 ) ) 

83  IF<J.GT.ISG(M))COTO  82 

IF(DATA(J).LT.XM(10*ID+M))THEN 
XM( 1Q*ID+M)-DATA( J) 

ENDIF 
J-J+2 
GOTO  83 

82  CONTINUE 

DO  44  J-lOdD-H.lO'ID+NSEG 
WRITE( 6 , 48 )  J 
DO  45  1-1.20 
DEEPd.J)-Cd.ID)-XM(J) 

IF ( DEEP (I . J ) . LT . 0 ) THEN 

DEEP (I ,  J)-0 

ENDIF 

WRITER  6. 25)  C(I,ID).AiT.J) .0(1. J) 
WRITE( 6 . * ) 


45  CONTINUE 

44  CONTINUE 

GOTO  49 

43  WRITE(6 , 24 )VS 

DO  46  1-1,20 
DEEP (I , ID)-C(I , ID)-ELO 
WRITE  (6.25)  C(I,ID),A(I.ID),Q(I.ID) 

WRITE  (6,*) 

46  CONTINUE 
RETURN 

C  COMPUTE  Z  FLOW  IN  EACH  SEGMENT 

49  DO  50  I-10*ID+1,10*ID+NSEG 
DO  50  J-1,20 

TQ( J, ID)-TQ( J, ID)+Q(J, I) 

50  CONTINUE 

DO  57  I-l.NSEG 
II«10*ID+I 
WRITE(6, 59)11 
WRITE{6, *) 

DO  57  J-1,20 

PERQ(J,II)-Q(J,II)/TQ(J.ID) 

IF(J.EO.l)  THEN 
PERQ( J , II )-0 
END  IF 

IF(PERQ(2. II ) .EQ. 1 . 0)THEN 

PERQ( 1,11 )-l . 0 

ENDIF 

IF ( ICODE .GT . 0 )THEN 

C(J,ID)-C(J,ID)*0.3048 

ENDIF 

WRITEC6, 60 )  C( J, ID) , PERQ( J, II ) 

WRITE(6. *) 

IF ( ICODE. GT.0) THEN 

C(J,ID)-C(J,ID)/0.3048 

ENDIF 

57  CONTINUE 

RETURN 

24  FORMAT (1H0.T4 2, 'RATING  CURVE  VALLEY  SECTION  ‘ ,F5. 1/T46, 'WATER' .Tj6. 
& ' FLOW' , T66 . 'FLOW1 /T45, ’SURFACE' .T56. 'AREA' ,T66, 'RATE' /T46. 'ELEV' . 
&T56 , ’SO  FT ' , T66 . ‘CFS’ ) 

48  FORMAT (1H0.T4 2. ‘RATING  CURVE  FOR  SEGMENT  ’, 15 . 1/T46 . ’WATER 756 , 

&’ FLOW' , T66 , ’ FLOW ’ / T4  5 , 'SURFACE' ,T56. ’AREA’ ,T66. ‘RATE’ /T46. ’ELEV’ . 
*T56, ’SQ  FT’ ,T66. ’CFS' ) 

31  FORMAT (1H0.T4 2. 'RATING  CURVE  VALLEY  SECTION ’. F5 . 1/T46 . 

&’ WATER’ ,T56, ’FLOW’ , T66 , ’FLOW’ /T45. ’SURFACE’ ,T56. ’AREA’ . 

&T66 , ’RATE’ /T4 6. 'ELEV' , T56. ’SQ  M’ ,T66. 'CMS' ) 

47  FORMAT(1HO.T4 2. ’RATING  CURVE  FOR  SEGMENT  ’.I5.1/T46. 

& ’ WATER ’ , T56 , ’FLOW’ .T66. *  FLOW’ /T45. 'SURFACE' ,T56, ’AREA’ , 

&T66 , ’RATE’ /T46 . 'ELEV' ,T56. 'SQ  M' ,T66. ’CMS’ ) 

25  FORMAT  ( 40X . F10 . 2 . 2F10 . 1 ) 

32  FORMAT  (40X.3F10.2) 

59  FORMAT  (1H0,T*7.’Z  DISCHARGE  IN  SEGMENT ’. 12 . 1 /T46 . 

&’ ELEV’ , T55 . 'PERCENT' ) 

60  FORMAT (40X.F10.2.2F10.3) 

70  FORMATdHO.T  10. 'ERROR  -  NEED  FLD  PLAIN  SEG  BOTH  S 
&IDES  OF  CHANNEL' ) 

170  FORMAT (1HO.T20. ’TURBULENT  EXCHANGE  METHOD ' . 2X . 1 5  1) 
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END 


SUBROUTINE  STT 

C  THIS  SUBROUTINE  STORES  A  DEPTH  -  FLOW  -  TRAVEL  TIME  TABLE. 

COfrWON / BLOCK 1 /  OCFS (300,6), DATA ( 3 1 1 ) , CFS ( 300 ) . CTBLE (50,11), 
&RAIN(300),ROIN(6), 

&A(20 , 70) ,Q(20 , 70) , DEEP (20 , 70 ) , ITBLE( 50 , 2) , DP(20 ) ,  SCFS( 20 ) , C(20 , 6) , 
&ZALPHA(20),IEND(6),DA(6) ,DIST(6) ,SEGN(6) ,DT(6) ,PEAK(6) ,ISG(6) , 

&NPU .NHD.NER, MAXNO , NCCtW ,  ICC , NCODE . TIME , KCODE ,  ICODE 
COWON / BLOCK2 /  PERQ(20 , 70 ) ,  TQ(20 . 6)  ,CC(20 )  ,LL(6) ,  INRC  ,LRC 

ID-DATA(l) 

REACH-DATA(2) 

MET1-DATAC5) 

IF (MET1 . EQ . 0 )GO  TO  2 
DATA (3 )-DATA( 3 ) /0 . 3048 
J-6 

DO  3  1-1,19 
DATA ( J ) “DATA (J)/0.3048 
DATA( J+l )-DATA( J+l ) /0 . 02832 
3  J-J+3 

2  XL-DATA(3) 

SLOPE-DATA ( 4 ) 

DIST ( ID  )-SLOPE*XL 
J-6 

DO  1  1-1,19 
DP(I)-DATA(J) 

SCFS(I )-DATA( J+l) 

CC(I)-DATA( J+2) 

*  J-J+3 

RETURN 
END 


SUBROUTINE  CMPTT 

C  THIS  SUBROUTINE  COMPUTES  THE  TRAVEL  TIME  AT  GIVEN 

C  DISCHARGE  RATES 

C 

C  IF  MULTIPLE  ROUTING  INVOKED,  COMPUTES  TRAVEL  TIME  TABLE  FOR 

C  THE  ONE  SEGMENT  SPECIFIED  *  OTHERWISE  ALL  SEGMENTS  TOGETHER 

C 

C  NOTE  FOR  MULTIPLE  RO"TINE  NEED  TO  REPEAT  THIS  ROUTINE  AND  ROUTE 
C  FOR  "EACH''  SEGMENT 

COWON / BLOCK  1  /  OCFS (300,6),  DATA(  3 1 1 ) .  CFS<  300 )  .  CTBLE  (50,11), 

&RAIN ( 300 ) . ROIN( 6 ) , 

&A(20. 70) ,0(20,70) . DEEP( 20 . 70 ) , ITBLE( 50,2).DP(20). SCFS( 20 ) ,C(20 . 6 ) . 
iZALPHA ( 20 ) , I END (6),DA(6),DIST(6>. SEGN(6) . DT( 6 ) , PEAK( 6 ) , ISG( 6 ) . 
fitNPU  .NHD.NER.  MAXNO .  NCOtW .  ICC  NCODE  .TIME.  KCODE .  ICODE 
COMMON /BL0CK2/  PERQ(20. 70 ) , TQ( 20 . 6) . CC( 20 ) . LL(6 ) . INRC.LRC 

ID-DATA(l) 
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REACH-DATA<2) 

NOVS-DATAO) 

IFOCCODE .  NE .  0 ) DATA (4  )«DATA( 4  )  /0 . 3048 
XL-DATA(4> 

SLOPE-DATA(5> 

DIST ( ID )-SLOPE*XL 
XLD36  -  XL  /  3600. 

MR-DATA<6) 

C  MULTIPLE  ROUTING 

INRC-DATA(7) 

C  RATING  CURVE  AT  TOP  OF  REACH 
LRC-DATAC8) 

C  RATING  CURVE  AT  DOWNSTREAM  END 

C  ZERO  ARRAYS 

DO  1  J-1,20 
DATA  (J)-0. 

1  CFS(J)-0. 

C  MULTIPLE  ROUTING  COMPUTATION 

IF(MR.LT. 1)GOTO  30 
WRITE(6,37) 

ID1-INRC 
GOTO  2 

30  ID 1-1 

C  FIND  RATING  CURVE  WITH  SMALLEST  MAXIMUM  FLOW  RATE 

2  QMIN-Q(20.ID1) 

MIN-ID1 

GO  TO  4 

31  ID1-LRC 
GOTO  32 

3  ID1-ID1+1 

32  IF  (QMIN-Q(20 , ID1 ) )  4.4.2 

4  IF(MR.LT.  DGOT033 
IFCIDl.EQ.INROGOTO  31 
IF ( ID1 . EQ . LRC )GOT05 

WRITE (6.*) ‘ERROR  only  two  r  curves  allowed  for  m. routing’ 
RETURN 

33  IF  (ID1-N0VS)  3.5.5 

5  I-I 
LL(ID)-0 

C  SET  SCFS  ARRAY  EQUAL  70  Q  ARRAY  OF  LOWEST  RATING  CURVE 
DO  6  J-2.20 
SCFS<I)0(J.MIN) 

IF (MR. LT . 1 )G0T0  6 
IF(PERQ( J.MIN) . LT . 0 . 001 ) THEN 
LL(ID)-LL(ID)+1 
END  IF 

6  r-r+i 

C  COMPUT  END  AREA  AND  DEPTH 

DO  9  ID1-1.NOVS 
IF(MR.LT.l)  GOTO  34 
IFCDl.EODTHEN 
ID1-INRC 
GOTO  34 
END  IF 
ID1-LRC 

34  M-l*LLcrD) 

N*2+LL( ID ) 

DO  36  >M.  19 
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DO  7  I-N.20 

IF  (Q(I , ID1)-SCFS< J) )  7,17,8 

7  CONTINUE 

17  DATA  <J)-A(I,ID1)+DATA(J) 

CFS(J)-DEEP(I.ID1)+CFS(J) 

GO  TO  36 

8  XY-(SCFS(J)-Q(I-1,ID1))/(Q(I,ID1>-Q(I-1,ID1)) 

DATA  ( J)*A!I~1 , ID1)+XY*(A(I , ID1)~A(I~1 , ID1) )+DATA( J) 

CFS( J)«DEEP(I-1 , ID1)+XY*(DEEP( I . ID1 )-DEEP(I-l, ID1) )+CFS{ J) 

36  CONTINUE 

IFCMR.LT.l)  GOTO  9 
IFCID1.EQ.LROGOTO  35 
ID1-1 

9  CONTINUE 

35  XNOVS-NOVS 

IF  C ICODE . EQ . 0 )GO  TO  19 
WRITE (6,20 )REACH 
GO  TO  21 

19  WRITE (6,13 )R£ACH 

ID1-MIN 

21  DO  10  I-M.19 

AVAREA  -  DATA  (I)  /  XNOVS 
DP  (I)  -  CFS(I)  /  XNOVS 
S  -  AVAREA  *  XLD36 
CC(I)-S/SCFS(I) 

IF ( SCFS ( I ) . EQ . 0 )  THEN 

CC(I)-0 

END  IF 

IF ( ICODE. EQ.O)GO  TO  24 
DP1-DP!I)*0.3048 
SCFS1-SCFS ( I ) *0 . 02832 
WRITE ( 6 , 14 )DP1 , SCFSl , CCC I ) 

GO  TO  10 

2U  WRITE ( 6 . 14 )DP! I ) , SCFS ( I ) ,CC( I ) 

10  CONTINUE 

C  PUNCH  CODE 

IFCNPU)12. 12,25 
25  IF! ICODE. EQ.01GO  TO  11 

XL1-XL*0 . 3048 

WRITE  ( 7 . 22  UD  .  REACH  ,  XL  L  .  SLOPE  .  ICODE 

DO  23  1-1, 19 

DP1-DP!  I )  *0 . 304  8 

SCFS1-SCFS! I )*0 . 02832 

WRITE! 7 , 26) DPI . SCFSl ,CCC I ) 

23  CONTINUE 

RETURN 

11  WRITE! 7, 1 5 ) ID . REACH , XL , SLOPE . ICODE 
WRITE  !  7  , 16  )  !DP(I).SCFS!I).CC(I). 1-1.19) 

12  RETURN 
C 

13  FORMAT ( 1H0 .746,  ' TRAVEL  TIME  TABLE ' /T54 . ' REACH' . F5  !//7*6. 'WATER' .7 
&56 , ' FLOW* , T65 , ' TRAVEL ’ /T*6 , 'DEPTH* ,T56,  RATE* .766. 'TIME'  T*6.  FEET 
&* , T56 , *CFS ' , T66 , 'HRS*  > 

14  FORMAT  (40X.F10.2.F10  0.F10.2) 

15  FORMAT! 'STORE  TRAVEL  TIME  ,  721 .  '  ID-  *  ,  II .  T29.  'REACH  NO-*,F5  1 ,  T  4  *  . 

&  '  LENGTH-’  ,  F9 . 0  .  ’  FTVT21.  'SLOPE-*  ,F8 .6.  ’FT/FT*  .  'CODE-MI  72 

5.1,  'DEPTH!  FT)  ’  .735.  ’FLCW(CFS)*  .749.  '  TIME  (HRS )  '  > 

27  FORMAT' 1 HO. 7^6. 'TRAVEL  TIME  TABLE ' /T5* .' REACH *. F6  1  T^  WATER*.! 


&56. 'FLOW' , T65 , 'TRAVEL' /T46 , ’DEPTH’ .156, 'RATE' ,T56, 'TIME' /T46, 
&''METER"  , T56 ,  ’CMS’  ,  T66 HRS’  ) 

22  FORMAT ( ' STORE  TRAVEL  TIME’ , T21 ID-\I1. T29, 'REACH  NO-’ ,F5. l.Ti*. 

&’ LENGTH-'  ,F9.0,  *  M'/T21,  'SLOPE-' ,  F8 . 6 .  ’M/M’  , "CODE-", I1/T2 
&1, 1  DEPTH (M)‘ ,  T35, 'FLOW(CMS)' ,T49, 'TIME (HRS)’ ) 

16  FORMAT  (T21.F7.2.F15.2.F15.3) 

26  FORMAT  (T21 ,F7 . 2 , 2F15 . 3 ) 

37  FORMAT (1HO.T24, ’MULTIPLE  ROUTING  INVOKED' ) 

END 


SUBROUTINE  ROUTE 

C  THIS  SUBROUTINE  ROUTES  A  HYDROGRAPH  THROUGH  A  REACH  WITH  THE 
C  NEW  VSC  METHOD  OF  FLOOD  ROUTING  THIS  METHOD  ACCOUNTS  FOR  THE 
C  VARIATION  IN  WATER  SURFACE  SLOPE. 

C 

C  IF  MULTIPLE  ROUTING  INVOKED  -  COMPUTES  PROPORTION  INFLOW 
C  FOR  ONE  SEGMENT 

C 

C  BUT  -  ONLY  ROUTES  ONES  SEGMENT  AT  A  TIME 

C  REPEAT  TRAVEL  TIME  TABLE  AND  ROUTE  COMMANDS  FOR  EACH  SEGMENT 
C  AND  ADD  OUTFLOWS 

COMMON / BLOCK 1 /  OCFS(3QO ,6) ,DATA(311) ,CFS(300) ,CTBLE(50 , 11) , 

&RAIN ( 300 ) , ROIN ( 6 ) ,  ' 

&A(20 . 70) , Q(20 , 70) , DEEP (20 , 70) , ITBLE(5G ,2) , DP (20) ,SCFS(2Q) ,C(20 , 6) . 
&ZAL PEA ( 2 0 ) . I END < 6 ) , D A ( 6 ) . D I S T < 6 ) , S EGN ( 6 > , DT ( 6 ) , PEAK ( 6 ) , I SG ( 6 ) , 

&NPU , NHD , NER . MAXNO . NCCfrM , ICC , NCODE .TIME, KCODE , ICODE 
COTWON/BLOCK2/  PERQ(20 , 70 ) ,TQ(20 , 6 ) ,CC(20> ,LL( 6) , INRC.LRC 
DIMENSION  DOCFS(300 . 6) 

C  New  variables  used 

C  DOCFS  copy  inflow  array  (preserve  original  OCFS) 

C  P  percentage  of  OCFS  (multiple  routing) 

C  DTDT  copy  inflow  time  increment  (preserves  DT(IDH)) 

C  IM  increments  DT 

ID-DATA(l) 

NHD-DATA(2) 

IDH-DATAO) 

DTUD)-DATA<4> 

DTDT-DT ( IDH ) 

DA(ID)-DA(IDH) 

M-IEND(IDH) 

MR-DATA(5) 

C  CHECK:  IF  M. ROUTING  ID.NE.IDH 

IF(MR.GT.O . AND. ID. EQ. IDH) THEN 

WRITE (6.*)' ERROR  -  FOR  M. ROUTING  ZD  MUST  NOT  BE  SAME  AS  IDH* 

RETURN 

ENDIF 

C  MULTIPLE  ROUTING  INCLUDED 

C  SET  UP  DUWff  ARRAY 

DO  51  1-1 , IEND(IDH) 

DOCFS ( I . IDH )  “OCFS ( I , IDH) 

MULTIPLE  ROUTING 

COMPUTE  DISTRIBUTED  FLOW  IN  SEGMENT 
IF (MR . LT . 1 )GOTO  50 


51 

C 

C 
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1 

II-INRC/10  ] 

NN-INRC-10*II  1 

IF  C ICODE .GT.O )GOTO  53 

WRITE (6 , 60 )NN  | 

GOTO  55 

53  WRITE ( 6 , 61 )NN  I 

55  TM-TIME-DTDT  j 

JJJJ-0  ^ 

DO  52  J-l.IEKD(IDH) 

TM-W-DTDT 

DO  56  K-2,20  < 

IF(DOCFS(J.IDH)-TQ(K,II))57, 58.56 

56  CONTINUE 

WRITE<6,*) 'FAILED  -  RATING  CURVE  EXCEEDED’ 

RETURN 

58  DOCFS{ J , IDH)*PERQ( J , INRC)*DOCFS( J , ID8) 

GOTO  5* 

57  ST-C(K.II)“(((TQ(K.II)-D0CFS(J.IDH))*(C(K.II)-C((K-1)III)))/(TQ 
&(K,II)-TQ<K-l.IT>M 

P-PERQ<K,INRC)-(((C(K,II)-ST)*(PERQ(K,INRC)-PERQ((K-1), 

&INRC)))/(C(K.II)-C((K-1).II))) 

DOCFS( J , IDH )-P*DOCFS ( J , IDH) 

5*  IF(DOCFS(  J ,  IDH)  .EQ.0)THEN 

IF ( J J J J ■ EQ . IEND ( IDH ) ) THEN 
WRITE(6, *)  'NO  FLOW  IN  SEGMENT’ 

RETURN 
END  IF 
GOTO  52 
END  IF 

IF ( ICODE. GT.O) THEN 

DOCFS ( J . IDH )-DOCFS (J . IDH ) *0 . 0283 168 

END  IF 

WRITE(6, 59)  TM , P , DOCFS ( J , IDH ) 

WRITE ( 6 , * ) 

IF (ICODE. GT.O) THEN 

DOCFS ( J . IDH ) -DOCFS ( J , IDH )/ 0 . 0283 168 

ENDIF 

52  CONTINUE 

C  IF  ID  AND  IDH  ARE  EQUAL.  ADD  1  TO  IDH 
50  IF(MR.LT. 1)THEN 
LL(ID)-0 
ENDIF 

IF  (ID-IDH)  3,1,3 

1  IDH-IDH+1 
DO  2  1*1, M 

2  DOCFS(I.IDH)«DOCFS(I,IOH-1) 

DT(IDH)-DTUDH-l) 

PEAK  IDH  )-PEAX(IDH- 1 ) 

3  NERRT-0 
PEAK (ID)  -  1. 

RO  -  0. 

N-19 

OCFS(1,ID)-0. 

S  -  0. 

T1  -  CC(1) 

J-l 

GUES  -  1. 
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A 
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CFS(l>-0. 

C  IF  ROUTING  INTERVAL  IS  NOT  EQUAL  TO  TIKE  INCREMENT  OF  INFLOW 
C  HYDROGRAPH ,  INTERPOLATE 

IF  <DT<ID)-DT(IDH>)  3.15,4 

4  TID-DT(ID) 

TIDH-0 . 

DO  7  1-2. M 
TIDB-TIDH+DT ( IDH ) 

IF  (TID-TIDB)  6.5.7 

5  J-J+l 
CFS(J)-DOCFSd.IDB) 
riD-TID+DTaC) 

GO  TO  7 

6  J-J+l 

CFS( J )-DOCFS(I-l, IDB)+( (TID-TIDH+DT ( IDH) )/DT( IDH ) )*  (DOC 
&FS(I, IDH)“DOCFS(I“l. IDH) ) 

TID-TID+DTCID) 

7  CONTINUE 
JO  TO  13 

8  TIDH-0. 

TID-DT(ID) 

DO  12  1*2. M 
TIDH-TIDH+DTdDH) 

9  IF  (TIDH-TID)  12.10.11 

10  J-J+l 

CFS ( J ) -DOCFS ( I . IDH ) 

TID-TID+DT( ID) 

IF  CJ-300)  12,13,13 

11  J-J+l 

CFSC J)-DOCFS< I- 1 . IDH)+( (TID-TIDH+DT ( IDH ) ) /DT ( IDH ) ) *(UOC 
&FS ( I , IDH ) -DOCFS ( I  * 1 . IDH ) ) 

TID-TID+DT(ID) 

IF  (J-300)  9,13,13 

12  CONTINUE 

13  DT(IDH)-DT(ID) 

W-J 

DO  14  1-2, M 

14  DOCFS(I, IDH)-CFS(I) 

C  IF  INFLOW  IS  ZERO.  SO  IS  OUTFLOW 

15  DO  16  L-2.M 

IF  (DOCFS (L, IDH) )  16.16.4g 

16  OCFS(L . ID )-0 . 

C  ROUTE 

49  DATA  (L-l )  -  0. 

DO  42  I-L.300 
IF  (I-M)  18,18,17 

17  DOCFS ( I , IDH ) -DOCFS (I -  1 . IDH ) * . 9 

18  AVIN- ( DOCFS ( I , IDH ) + DOCFS ( I - 1 . IDH ) ) / 2 . 

SIA  -  AVIN  ♦  S 

J-l 

C  DETERMINE  DEPTH  AND  TRAVFr,  TImt  of  INFLOW 

IF  ( DOCFS ( I . IDH) -SCFS( 1+LL( ID) ) )  19.23.20 

19  D 12  -  (DOCFS(l , IDH)  /  SCFS( 1+LLIID) ) )  *  DP(1*LL(ID>) 

TI2  -  CC(  1+LLdD) ) 

GO  TO  25 
JJJ-2 

IF(LL( ID) .GT . 0 ) THEN 
JJJ-LL(ID)+2 
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END  IF 

DO  21  J-JJJ.N 

IF  (DOCFS(I , IDH)-SCFS( J) )  24,23,21 

21  CONTINUE 

IF  (NERRT)  22,22,36 

22  WRITE  (6,46) 

NERRT- 1 

GO  TO  36 

23  DI2-DP<J) 

TI2  -  CC(J) 

GO  TO  25 

24  RATICH(DOCFS(I.ID0)-SCFS(J*1))/(SCFS(J)-SCFS(J>1)) 
DI2-DP(J-l)+RATIO*(DP(J)-DP<J-l)) 
TI2-CC(J-l)+RATIO*(CC(J)-CC(J-l)) 

25  DO  35  IT-1.10 
J-l 

C  DETERMINE  DEPTH  AND  TRAVEL  TIME  OF  OUTFLOW 

IF  (GUES-SCFS(ULLUD)))  26,29.27 

26  D02  -  (GUES  /  SCFS( 1+LL(ID) ) )*  DP( 1+LL(ID> ) 

T02  -  CC(1+LL(ID) > 

GO  TO  31 

27  DO  28  J-JJJ.N 

IF  (GUES-SCFS(J))  30,29,28 

28  CONTINUE 
J-N 

29  D02-DP(J) 

T02-CC(J) 

GO  TO  31 

30  RATIO- ( GUES -SCFS( J-l ) ) / (SCFS( J)-SCFS( J-l ) ) 

D02-DP( J-l)+RATIO*(DP( J)-DP( J-l) ) 

T02-CC ( J- 1 ) * RATIO* (CC(J)-CC(J-l)) 

C  FIND  WATER  SURFACE  SLOPE 

3 1  DDD-DIST ( ID ) / (DIST(ID  J+DI2-D02 ) 

IF  (DDD-.01)  32.32,33 

32  GUES-DOCFS(I-l.IDH) 

GO  TO  35 

33  T2  -  .5  *  (TI2  ♦  T02) 

T2-T2*SQRT(DDD) 

T  •  T1  ♦  T2 

C  COMPUTE  ROUTING  COEFFICIENT 

COEF  -<2.  *  DT(ID) )  /  (T+DT(ID) ) 

02  -  COEF  *  SIA 
TRY1  -  GUES 
RATIO-02/ (GUES+ . IE-20 ) 

DIFF-ABS(1. -RATIO) 

C  TEST  FOR  CONVERGENCE 

IF  (DIFF-0.001)  37,37.34 

34  GUESH32 

35  CONTINUE 
OCFS(I.ID)-DATA(I-l)*SIA 
DATA(I)  -  DATA(I-l) 

WRITE  (6.47)  I.OCFSU.ID) 

GC  TO  3B 

36  OCFS(I.ID)-DATA(I-l)*SIA 
DATA(I)  -  DATA(I-l) 

GO  TO  38 

37  OCFS(I.ID)-02 
DATA  (I)  -  COEF 
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C  COMPUTE  MEW  STORAGE 

38  S  -  SIA  -  OCFSU.ID) 

T1  -  T2 

RO  -  RO  +  OCFS  Cl, ID) 

IF  (OCFS (I. ID)  -  OCFS(I-l.ID))  39,40.40 

39  IF(OCFSCI.ID)  -1.)  43,43,42 

40  IF(OCFS(I, ID)  -  PEAK(ID))  42,42,41 

41  PEAX(ID)-OCFS  (I, ID) 

42  CONTINUE 
1-300 

43  IEND(ID)-I 

ROINC ID )  -  RO*DT(ID)*3600 
C  PUNCH  CODE 

IF  (MPU)  45.4S.44 

44  WRITE  (7,48)  ID.NHD, IDB.DT(ID) 

45  DT(IDH)-DTDT 
RETURN 

C 

46  FORMAT (1H0 ,  ’TRAVEL  TIME  TABLE  EXCEEDED ’ ) 

47  FORMAT (T 10, 'PROBLEM  FAILED  TO  CONVERGE  AFTER10  ITERATIONS .  CONVEKG 
&ENCE  WAS  FORCED. */T20, 'OUTFLOW  NUMBER  -  '.14,’RATE  «',F10.2) 

48  FORMAT C  'ROUTE' ,T21, ' ID-’ , I1.T29. 'HYD  NO-’ , I3.T45, ’ INFLOW  ID-* , I 
&1.T65, ’DT-' , F8 . 6, ’HRS') 

60  FORMAT (1H0.T40. 'INFLOW  FOR  SEGMENT 15 . 1/T30 HOURS T40 , 

&' PERCENT’ ,T52, 'CFS') 

61  FORMAT (1H0.T40. 'INFLOW  FOR  SEGMENT ’, 15 . 1/T30 HOURS ’ , 

&T40, ’PERCENT’ ,T52, 'CUMECS’ ) 

59  FORMAT ( 25X ,  F10 . 3 , 2F10 . 3 , 3FK  3) 

END 


SUBROUTINE  RES VO 

C  THIS  SUBROUTINE  ROUTES  A  HYDROGRAPH  THROUGH  A  RESERVOIR  WITH  THE 
C  STORAGE- INDICATION  METHOD. 

COWON / BLOCK  1  /  OCFS (300, 6)  ,DATA(311) , CFS (300)  .CTBLEC50, 11), 
&RAINC300 ) , ROINC 6) , 

&AC20.70) ,0(20, 70),DEEP(20. 70). ITBLEC 50, 2) ,DP(20) ,SCFS(20) , 0(20.6). 
&ZALPHAC20) , IEND(6) ,DA(6) ,DIST(6) ,SECN(6) ,DT(6) , PEAK(6) , ISGC6) , 

&NPU .NHD.NER, MAXNO . NCOW , ICC . NCOPE . TIME , KCODE , ICODE 
CCWON / BL0CK2 /  PERQ(20, 70) ,  TQ'2«. ,  6)  ,CC(20)  ,LL(6) ,  INRC,  LRC 

ID-DATA(l) 

NHD-DATA ( 2 ) 

IDH-DATAO) 

NERES-0 

DT( ID)-DT(IDH) 

RO  -  0. 

DA( ID)-DA(IDH) 

PEAK(ID)  -  1. 

J-l 

1-4 

C  REMAINING  DATA  ARE  FLOW  AND  STORAGE  VALUES 
IF (KCODE. EQ.0)GO  TO  25 
DATA(I)-DATA(I)/0. 02832 
DATA (1+1 )«DATA ( 1+ l ) / l . 2 1 968 


25  SCFS ( J ) "DATA ( I ) 

STOR£l-DATAd+l)*12.2 
STORE-STORE 1 

C  COMPUTE  STORAGE  COEFFICIENT  ARRAT  c 

1  CC(J)-(SCFS(J)/2. )+(STOR£/DT(ID) ) 

I-I+2 

J-J+l 

IF  (J-20)  2,2,3 

2  IF  (KCODE ,  EQ .  0  )GO  TO  26 
DATA<I)-DATAd)/0. 02832 
DATAd+ll-OATAd+D/l. 21968 

26  SCFS(J)-DATAd) 

STORE-DATA! 1+1 )* 12.1 
IF  (SCFS(3>“ .001)  3,3,1 

3  H-J-l 
OCFSd.  ID)-o. 

S-STOR£l/DT(ID) 

C  ROUTE 

00  13  1*2,150 
IF  (I-IEND(IDH))  5.5,4 
*  OCFSd. IDH1-0.0 

5  AVIN-(0CFSd.IDH)4OCFSd-l.IDH>  >/2. 

SIA-S+AVIN 

c  DETERMINE  PROPER  C 

DO  6  J-l.N 

IF  (SIA-CC(J))  10,9.6 

6  CONTINUE 

IF  (NERES)  ),2,a 
I  WRITE  (6,19) 

NERES- I 

8  RESC-SCFS  <  N ) /CC  <  N ) 

C  COMPUT  OUTFLOW 

OCFSd. ID)-RESC*SIA 
GO  TO  11 

9  OCFSd,  ID)-SCFS(J) 

GO  TO  11 

10  OCFSd, ID)-SCFS(J-1)TUSIA-CC(J-1))/(CC(J)-CC(J 

s  -1)))*(SCFS(J)-SCFS(J-1)) 
c  DETERMINE  NEW  STORAGE 

11  S-SIA-OCFSCI , ID) 

RO  -  RO  *■  OCFSd,  ID) 

IF  (OCFSd,  ID)-0CFS(I-1.  ID))  12,13.13 

12  IF  (OCFSd, ID)-1.  )  16,16,15 

12  IF(OCFSd.ID)  -  PEAKdD))  15,13.14 

11  PEAKdD)  -  OCFSd.  ID) 

15  CONTINUE 
1-150 

16  IEND(ID)-I 

ROIN(ID)  -  RO  •  DT(ID)»3600 
C  PUNCH  CODE 

IF  (NPU)  18.18,17 

12  II-2-N+3 

IF ( ICODE . EQ , 0 )GO  TO  22 
WRITE(7,24)ID.NHD,  IDH. KCODE 
DO  23  1-5,11,2 
DATA(I)-DATAd)  *0.02832 
DATAd  +  l)-DATAd*l)*  1.21968 

23  CONTINUE 
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WRITE ( 7 , 2?) (DATA< I ) , 1*5, II ) 

RETURN 

22  WRITE( 7 , 20 ) ID , NHD , IDH , ICODE 
WRITE  (7,21)  (DATA( I ) , I"5, II ) 

18  RETURN 
C 

19  FORMAT  (1H0, 3 3HSTORAGL -DISCHARGE  TABLE  EXCEEDED.) 

20  FORMAT  (  ’ROUTE  RESERVOIR  T21 ,’  ID-', II ,  T29,  '  HYD  NO-',I3 ,  T*2  IMF 

&LOW  ID"  ’ , I 1 , T60 , "CODE"" , I 1  /T21 . ’OUTFLOW(CFS ) ' , T37 , ' STOR 

&AGECAC  FT) * ) 

2k  FORMAT!  'ROUTE  RESERVOIR' ,T21. 'ID-’ ,11, T29, ’HYD  NO-\I3,T42, 'INF 

&LOW  ID" * ,11, T60 , "CODE-" , I 1  /T21, ' OUTFLOW ( CMS ) ' ,T37, 'STOR 

&AGEC 1000CU  M)’ ) 

21  FORMAT  (T21 ,F10. 1.F13. 1) 

27  FORMAT  (T21 , F10 . 2 , F13 . 2) 

END 


SUBROUTINE  ERROR 


C  This  subroutine  determines  the  error  standard  deviation  and  the  peak  flow 
C  error  for  2  hydrographs  (original  program  retained). 

C  Assumes  that  measured  is  ID1 

C  In  addition,  10  other  measures  of  goodness  of  fit  are  calculated. 

C  All  Indicies  are  printed  out  in  metric  units. 

C0M«N  /  BLOCK  1  /  OCFS  (300,6),  DATA  (  3 1 1 ) .  CFS  (  3 00  ) ,  CTBLE (50,11). 
&RAIN(300),ROIN(6), 

&A(20,70),Q(20,70), DEEP (20 , 70 ) , ITBLE( 50 , 2) , DP(20) , SCFS(20 ) , C{ 20 . 6 ) , 
&ZALPHAC  20 ) , IEND(6 ) ,DA(6) , DISTC6) , SEGN(6 ) , DT(6 ) , PEAK(6 ) , ISG(6 ) . 

&NPU . NHD , NER , MAXNO , NCCW1 , ICC , NCODE , TIME , KCODE , ICODE 
ID1"DATA( 1) 

ID2-DATAC2) 

SSE-0 . 

WRITE(6 . 21 ) 

21  FORMAT (1H0.T33, 'TIME' ,T55, 'FLOW  l',T?6, 

&  'FLCW  2 ’ , T95 , ' ERROR’ /T34 , 

&  ’HRS’ ,T57. ’CMS’ ,T78, ’CMS’ ,T97, 'CMS' ) 

22  J-l 

C  If  the  time  Increments  are  not  equal,  interpolate. 


IF  (DT(IDl)-DT( ID2) )  1,8,2 

1  L-ID1 
K-ID2 
GO  TO  3 

2  L-ID2 
K-ID1 

3  M-IEND(L) 

TID-DT(K) 

TIDH-0. 

DO  6  1*2 ,M 
TIDH-TIDH+DT ( L ) 

IF  (TID-TIDR)  5,4,6 
k  J-JM 

CFS  ( J  )"OCFS  ( I ,  L ) 
riD"TID+DT(K) 


297 


GO  TO  6 

5  J-J+l 

CFS(J)MXFS(I-1.L)  +  (<TID-TIDB+DT(L))/DT(L)>*<0CFSU,L)-0CFS(I-1.L> 
&> 

TID-TID+DT(K) 

6  CONTINUE 
IEND(L)«J 
DT(L)-DT(K) 

DO  7  1-2. J 

7  OCFS(I.L)-CFS{I> 

8  IF  ( IEND( ID1)”IEND(ID2) )  9,9,10 

9  M-IEND(IDl) 

GO  TO  11 

10  M-IEND(ID2) 

11  T2-TIME 

IF  ( KCODE . EQ . 0 ) THEN 
DO  997  I-l.M 

OCFS ( I . ID1 ) -OCFS ( I , ID1 ) * . 02832 
997  OCFS(I , ID2 )-OCFS(I , ID2)* . 02832 

END  IF 


C  Determine  error  -  original  method 
DO  12  I-l.M 

ERR-OCFS { I . ID 1 ) -OCFS ( I , ID2 ) 

WRITE(6 , 16)T2 ,OCFS( I , IDl) ,OCFS( I , ID2 ) , ERR 
l6  FORMAT  (6X.F12.3.3F12.0) 

25  T2-T2+DT(ID1) 


C  Sum  of  squares  of  error 

12  SSE-SSE+ERR*ERR 
XM-M 


C  Error  variance 
EVAR-SSE/XM 


C  Error  standard  deviation 


ESDEV-SQRT(EVAR) 


C  Percent  error  for  peak  discharge 


ERPK-ABS ( PEAK ( I D 1 ) -  PEAK ( I D2 ) ) 

FCTER-(ERPK/PEAK( IDl ) )*100. 

C  Other  goodness  of  fit  calculations... 


SUMO 1-0. 
SUMO-O . 
SUM 1-0. 
SUM2-0 . 
SUM3-0 . 
SUMA-0 . 
SUM5-0 . 
SUM6-0 . 
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SUM7-0. 

SUM8-0. 

SUM9-0 . 

SUM1 0-0. 

SUMU-0 . 

SUM 12-0 . 

DO  77  I-l.M 

ERR-OCFS ( I . ID1 ) -OCFS ( I , ID2 ) 

IT (OCFS (l , 101 ) . EQ. 0 . 0 .AND. OCFS (I , ID2 ) . NE. D . 0)THEN 
LOGERR-ALOG ( OCFS ( I , ID2 ) ) 

ELSE  IF(OCFS(I,ID1>.NE.O.O.AND.OCFSU.ID2).EQ.O.O)THEN 
LOGERR-ALOG (OCFS (I, ID1)) 

ELSE  IF(0CFS(I,ID1)  .EQ.O. 0. AND. OCFSC I. ID2)  .EQ.O. OTHER 
LOG  ERR -0 . 

ELSE 

LOGERR-ALOG ( OCFS (1 , ID 1 ) ) -ALOG ( OCFS (I , ID2) ) 

ENDIF 

sumo -ocfs  ( r ,  ro  1 )  ■*  sumo 

SUMO 1-OCFS (I . ID2 ) +SUM0 1 

SUM1-ERR+SUM1 

SUM2-ERR*#2+SUM2 

SUM3-LOGERR**2+SUM3 

IF (OCFS (I , ID1) .EQ. 0 • >OCFS(I . ID1)-1 . 

SUM4-( (ERR/OCFS(I . ID1 ) )**2>+SUM4 
77  CONTINUE 

DO  13  1-2, M 

DIFF1-OCFS(I,ID1)-OCFS(I-1,ID1) 

DIFF2-OCFS (I . 1D2 ) -OCFS (1-1,102) 

SUM5-( (DIFF1-DIFF? )**2)+SUM5 
SUM7-DIFF1+SUM7 
IF(DIFF1.EQ.0. )DIFF1-1. 

SUM6- ( ( (DIFF1“D1FF2)/DIFF1)#*2)+SUM6 
13  CONTINUE 


SltMuAN-SUMOl/M 

OBSMEAN-SUMO/M 

DIFFM1-SUM7/M 

DO  1*  1-2. M 

SUMS- ( ( (OCFS ( I . ID 1 ) -OCFS ( I - 1 , ID 1 ) ) -DIFFM1 ) **2 ) +SUM8 
SUM9- ( ( ( (OCFS ( I , ID1 ) -OCFS ( r- 1 . ID 1 ) ) /DIFFM1 ) - 1 ) **2 ) +SUM9 

H  CONTINUE 

DO  73  I-l.M 

SUM1 0- ( ( OCFS ( I , ID 1 ) - OBSMEAN ) **2 >+SUMl 0 
SUM1 1- ( ( ( OCFS ( I . ID 1 ) /OBSMEAN ) - 1 ) **2 ) +SUM1 1 
SUM12-( (OCFStl . ID2)-SI>WEAN)**2)+SUM12 
73  CONTINUE 

SDM-SQRT ( SUM1 0 / (M- 1 ) ) 

SDS-SQRT ( SUM1 2 / ( M- 1 > ) 

DO  115  1-1 ,M 

115  SUMl3-( (OCFS (I , ID1 ) -OBSMEAN) /SDM) *( (OCFS (I , ID2 ) - 

&  SIhM£AN)/SDS)*SUMl3 
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i 

I 


f 

1 

( 

i 

t 

t 


A 


OF 1-SUM 1 
OF2-SUM2 
OF3-SUM3 
OF 4 -SUM* 

OF5-SUM5 

0F6-SUM6 

OF7-SUM2/SUM10 

OF8-SUM4/SUM11 

OF9-SUM5/SUM8 

OF10-SUM6/SUM9 

AM-M 

0F11-(1.  /AM;*SUZil3 
WRITEC6 , 95) 

95  FORMAT ( 1H0, 10X, ' . . 

WRITE(6, 50) 

50  FORMAT ( 1 5X , *  MEASURES  OF  FIT  ’//) 
WRITEC6, 91) 


91 

FORMAT (1 OX, . . 

WRITE(6 , 51 )0F1 

51 

FORMAT (1 OX. ’SUM  OF  ERRORS 

WRITE(6, 52)OF2 

’ . F20 . 5) 

52 

FORMAT (1 OX. 'OLSO 

WRITEC6, 53)OF3 

’ , F20  5) 

53 

FORMAT (1 OX. ’LOG  LSQ 

WRITE (6 , 54 )OF* 

’ . F20  5) 

54 

FORMATC1 OX, 'RELATIVE  ERROR 

WRITE{6, 55)OF5 

’ .F20.5) 

55 

FORMAT ( 10X . * ABS  ERROR  -  DIFF 

WRITE ( 6 , 56 )OF6 

’ , F20 . 5) 

56 

FORMAT ( 10X . ’ REL  ERROR  -  DIFF 

WRITE(6, 57 )OF7 

' .F20.5) 

57 

FORMAT (1 OX. ’ABS  ERROR/VAR 

WRITE (6,58 )OF8 

’  .  F2;  .  5 ) 

50 

FORMATUOX.  'REL  ERROR/VAR 

WRITE(6,59)OF9 

’ , F20 . 5 ) 

59 

FORMATdOX, 'ABS  ERRORCdiff )/VAR 

WRITE (6. 60 )OF 10 

’ . F20 . 5) 

60 

FORMATdOX. ’REL  ERROR(diff  )/VAR 

WRITE( 6 . 61 )OFll 

’ .F20. 5) 

61 

FORMAT ( 10X . ’ PEARSONS  r 

WRITE(6 , 92 )ESDEV 

’ ,F20  5) 

92 

FORMAT (10X . ’ ERR  STANDARD  DEV 

‘.T1ITE(6,93)PCTER 

’ . F20 . 5) 

93 

FORMAT (10X . ’ PEAK  Q  ERROR 

WRITE(6, 96) 

’ .F20.5) 

96 

FORMATdOX.  ’ . — . 

WRITE  (6.98) 

98  FORMAT  ( / / 10X  ’NOTE:  All  indicies  are  in  metric  units’) 

IF  (KCODE . EQ . 0 ) THEN 
DO  9969  1-1 ,M 

OCFS(I,ID1)-OCFS(I.ID1)/.02832 
9969  OCFS(I.ID2)-OCFS(I. ID2)/. 02832 

ENDIF 
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RETURN 

C 

END 


SUBROUTINE  SEDT 

C  THIS  SUBROUTINE  COMPUTES  THE  SEDIMENT  YIELD  FOR  A  FLOOD 

cotton/ BLOCKl /  OCFS (300.6). DATA( 3 H ) . CFS (300 ) . CTBLE (50.11), 
&RAINOOO) , ROIN(6) , 

&A( 20 . 70) ,Q(20 . 70 ) . DEEP( 20 . 70 ) . ITBLE( 50 . 2 ) . DP( 20 ) . SCFS (20 >,0(20.6), 
&ZALPHA(20) . IEND(6 ) , DA(6 ) , DIST(6 ) ,SEGN(6 ) ,DT(6 ) , P£AK(6 ) , ISG(6) , 

&NPU .  NHD ,  NER .  MAXNO .  NCOW .  ICC.  NCODE .  TIME .  KCODE .  ICODE 

ID-DATA ( 1 ) 

SOIL-DATA(2) 

CROP-DATA(3) 

CP-DATA(4) 

SL-DATA(5) 

WRITE(6, *)  ***  CHECK  THIS  IS  CORRECT  AREA’ , DA( ID) 

WRITE ( 6, •)• ESPECIALLY  IF  MULTIPLE  ROUTING  UTILIZED’ 

C  COMPUTE  SEDIMENT  YIELD 

X“ROIN( ID) #DA( ID )*53 . 333 * PEAK ( ID) 

S ED-95 . *X** . 56*SOIL*CROP*CP*SL 
IF (ICODE  EQ. 0 )GO  TO  5 
SED1-SED*0 . 9072 
WRITE(6 . 6 )SED1 
GO  TO  7 

5  WRITE  (6.3)  SED 

C  PUNCH  CODE 

7  IF(NPU )2 . 2 . 1 

1  WRITE  (7.4)  ID. SOIL. CROP. CP. SL 

2  RETURN 

3  FORMAT  (I0X.  ’SEDIMENT  YIELD  -  ’.  F10.1.  '  TONS’) 

4  FORMAT (  ' SEDIMENT  YIELD’ . T21 , ' ID- ’.II. T29 , ’ SOIL-’ . F5 . 3 ,T42 , ’CROP 

.F5.3.T57. 'CP-’ .F5.3.T70. 'LS-*  F5.3) 

6  FORMAT (10X. "SEDIMENT  YIELD-" . F 10 . 1 .  'METRIC  TON") 

END 


BLOCK  DATA 

C  BLOCK  DATA  SUBPROGRAM  UZED  TO  INITIALIZE  ZALPHA . CTBLE , IT BLE 
C  AND  NCC*M . 

COMCN  / BLOCKl  /  OCFS  (300.6),  DATA  C  3 1 1 )  .CFS(  300  5 .  CTBLE{  50,11). 

&RA IN( 300 ) ,ROIN(6) , 

&A(20 , 70). 0(20, 70). DEEP (20 . 70 ) , ITBL£( 50 , 2) , DP( 20 ) , SCFS ( 20 ) . C ( 20 . 6 ) . 
4ZALPHA ( 20 ) . I END (6). DA (6). DI ST (6). SEGN(6 ) , DT( 6 ) , PEAK( 6 ) , ISG( 6 ) . 

&NPU .  NHD  .  NER . MAXNO .  NCOM .  ICC .  NCODE .  TIME .  KCODE  .  ICODE 

COftON / BLOCK2 /  PER0( 20 . 70 ) . TQ(20 . 6) .CCl 20 ) . LL( 6 ) . INRC , LRC 


DATA  ZALPRA/THl. 1H2. IH3. 1H* . 1H5. 1H6. 1H7. IR8 . 1H9. 2H0, 2H 
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.18. , 1H-.1H  ,  1H  .1H  , 1H  , 1H  / 

DATA  NCOW/17/ 

DATA  CTBLE/ IBS , IBS , 1HR.1HC, IBP, IBP. 1HP . 18A , 1HS , 1BC , 1HS , 1HC , 1HR , 
&1HR , 1HE , IBS , 2HF, 33*1H  . 

&1HT  .  1BT ,  1SE  ,  1H0 . 1HR,  1HU,  1HL ,  1HD,  1HT  ,  1H0 . 1BT ,  1B0 . 1H0 , 1HQ . 1HR , 
&1BE,1BI.33*1B  , 

&2HAR , 2H0R . 2HCA . 2RMP . 2HIN , 2SSC , 2B0T , 2 HD  , 2H0R , 2HMP . 2H0R . 2HMP , 
&2HUT , 2HUT , 2HR0, 2HDI , 2HNI . 33*2H  , 

&2BT  , 2HE  , 2HLL , 2HUT, 2ST  , 2HH  . 2H  H.2HHY.2HE  ,2Bt’T,2H£  .2BUT, 
&2BE  ,2HE  , 2 SR  , 28ME , 2HSB , 33*2H  , 

&2B  , 2HHY, 2fl  H, 2 HE  , 2HHY , 2HBY , 2BYD , 2BD  , 2HRA . 2BE  , 2HTR , 2HE  . 

&2B  .2HRE.2HAN, 2HNT, 2H  ,33*2B  , 

&2H  , 2HD  , 2SYD , 2HHY , 2HD  , 2BD  , 2H  ,2H  , 2HTI . 2HRA , 2HAV , 2HTR , 

&29  , 2BSE , 2BAL , 2H  Y.2B  . 33*2H  , 

&2B  ,2H  . 2B  ,2HD  ,2H  , 2H  ,28  , 2B  , 2HNG , 2HT I , 2HEL , 2HAV , 

&2B  .2HRV.2BYS.2HIE.2H  .33*2H  . 

&8-2R  . 2H  C.2HNG.2H  T.2HEL.2H  , 2H0I . 2HIS , 2HLD , 34*2H  , 

&8-2H  .2HUR.2H  C.2HIM.2H  T.2H  , 2HR  ,36*2H  . 

68- 2B  , 2BVE , 2HUR , 2HE  ,2HIM,38*2B  , 

69- 2H  .2HVE.2H  , 2BE  .38*2H  / 

DATA  ITBLE/1, 2. 3. 4. 5. 6, 7.8.9.10.11.12,13, 14  .15.  *6. 17. 33*lH  . 

&4,  3 10. 310, 3 10. 4. 1.2. 4. 100, 310. 100, 8. 7, 25. 2.5.0,33MB  / 

EHD 


> 


*  EXAMPLE  datal  DATASET 

START  00.00  0  0  1 

*  RUNOFF  SUBCATCHMENT  <<06,  MARBACH 

*  SUBCATCHMENT  406 

STORE  HYD  ID-2  HYD  NO-406  DT-2.0  HRS  DA-57.1  SQ  MI 

BASEFLOW  -  154  CFS 

FLOW  RATES  -717  653  635  632  580  613 
907  1095  950  483  295  322  296  173  99  65  63 
67  43  27  18  13  10  7  5  4  3  2  85  185  98  42  18 
13  9  20  33  18  9  103  274  239  185  315  356  186 
141  168  99  53  31  23  101  197  107  147  238  131 
61  30  22  16  12  9  7  5  4  3  2  1  1 


•  ROUTE  RUNOFF  TO  HERMANNS PI EGAL 

*  COMPUTE  RATING  CURVE  FOR  MARBACH 

COMPUTE  RATING  CURVE  ID-1  IT-1  MR-1  VS-6  NO  SEGS-3 

MIN  ELEV-831.0FT  MAX  ELEV-869.4FT 
CH  SLP-0.006  FLDPN  SLP-0.0075 


N-0.05 

DIST-111 . 

6FT 

N—  .03 

DIST-14 1 . 

1FT 

N-0.05 

DIST-344. 

8FT 

DIST 

ELEV 

DIST 

ELEV 

FT 

FT 

FT 

FT 

0.0 

869.4 

32.8 

856.3 

40.0 

855.3 

52.5 

854.3 

65.6 

047.8 

80.4 

845. 

91.9 

841.9 

105.0 

837.9 

111.6 

836.9 

114.8 

835.3 

121.4 

831.0 

134.5 

831.0 

141.1 

836.6 

144.4 

836.9 

157.5 

839.6 

170.6 

840  2 

180.4 

845.5 

344.8 

859.6 

*  COMPUTE  RATING  CURVE  FOR  HERMANNS PI EGAL 
COMPUTE  RATING  CURVE  ID-2  IT-1  MR-1 

VS- 7  NO  SEGS-3  MIN  ELEV-687.7FT 
MAX  ELEV-702 . 1ft 


CH  SLP 

-  0.006 

FLDPL 

SLP-0.0075 

N-0.05 

DIST-118 

.4FT 

N— .03 

DIST-171 

.  6FT 

N-0.05 

DIST-220 

.  1FT 

DIST 

ELEV 

DIST 

ELEV 

FT 

FT 

FT 

FT 

0.0 

702.1 

49.  Z 

697.5 

105.6 

696.5 

118.4 

696.2 

124.7 

693.9 

132.2 

688.3 

133.5 

688.3 

136.2 

687.7 

147.6 

688.0 

154.2 

688.3 

.5.  5 

689.  0 

257.5 

689.6 

167.3 

695.5 

171.6 

696.2 

187.0 

696.2 

200.1 

696.5 

212.3 

696.2 

216.5 

698.5 

220.1 

722.1 

•  LEFT  SEGMENT 

COMPUTE  TRAVEL  TIME  ID- 3  REACH  NO- 4  NO  VS-2 

L-110564.3  FT  SLP-0.007  MR-l 
INRC-11  LRC-21 
ID-3  HYD  NO-141  INFLOW  ID-2 
DT-0 . 5  KRS  MR-1 


ROUTE 
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*  CHANNEL  SEGMENT 


COMPUTE  TRAVEL  TIME  ID-4  REACH  NO-4  NO  VS-2 

L-I10564 . 3  FT  SLP-0.007  KR-1 
INRC-12  LRC-22 

ROUTE  ID-4  HYD  NO-142  INFLOW-2 

DT-0.5  HRS  MR-1 

*  RIGHT  SEGMENT 

COMPUTE  TRAVEL  TIME  ID-5  REACH  N04  NO  VS-2 

L-I10564 . 3  FT  SLP-0.007  MR-1 
INRC-13  LRC-23 

ROUTE  ID- 5  HYD  NO-143  INFLOW-2 

DT-0.5  HRS  MR-1 

*  ADD  SEGMENTS 

ADD  BYD  ID-2  HYD  NO407  IDI-3  IDII-4 

ADD  HYD  ID-2  HYD  NO407  IDI-2  IDII-5 

*  _ _ _ __ _ _ _ ___ _ 

*  RUNOFF  SUBCATCHMENT  407 

*  SUBCATCHMENT  407 


STORE  HYD  ID-3  HYD  NO407  DT-2.0  HRS  DA-105.7  SO  MI 

BASEFLOW  -  83  CES 


FLOW  RATES-  1022  1155  1318  1403  1401  1320 
1380  1698  2019  1807  1261  840  775  717  550  389 
281  210  160  133  112  93  78  66  55  46  38  32  56 

187  248  197  122  75  47  53  155  203  161  124  217 

354  332  *58  618  602  448  41*  478  433  305  199 

143  154  157  191  456  576  *52  279  173  109  71  58 

48  40  14  28  24  '20  16  »  11  9  8  6  5  *  t  1  !  2  1 


*  OUTFLOW  HERMANNSPIEGAL 

m  DTD  ID-2  HYD  NO-407  IDI-3  IDII-2 

PRINT  HYD  ID-2  NPK-1  rOR-O  IN-0 

FINISH 


